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Preface 



A volume of this nature containing a collection of papers has been 
brought out to honour a gentleman - a friend and a colleague - whose work has, 
to a large extent, advanced and popularized the use of stochastic point 
processes. 

Professor Srinivasan celebrated his sixty first birth day on December 
16,1990 and will be retiring as Professor of Applied Mathematics from the 
Indian Institute of Technology, Madras on June 30,1991. In view of his 
outstanding contributions to the theory and applications of stochastic 
processes over a time span of thirty years, it seemed appropriate not to let 
his birth day and retirement pass unnoticed. A symposium in his honour and the 
publication of the proceedings appeared to us to be the most natural and 
suitable way to mark the occasion. The Indian Society for Probability and 
Statistics volunteered to organize the Symposium as part of their XII Annual 
conference in Bombay. We requested a number of long-time friends, colleagues 
and former students of Professor Srinivasan to contribute a paper preferably 
in the area of stochastic processes and their applications. The positive 
response and the enthusiastic cooperation of these distinguished scientists 
have resulted in the present collection. 

The contributions to this volume are divided into four parts: Stochastic 
Theory (2 articles). Physics (6 articles), Biology (4 articles) and Operations 
Research (12 articles). In addition the keynote address delivered by 
Professor Srinivasan in the Symposium is also included. 

A large number of individuals have helped to make this volume a worthy 
scientific monument for Professor Srinivasan. First among these are the 
authors and the referees. Several authors have added personal comments 
expressing their high esteem and respect for Professor Srinivasan as a 
researcher, colleague and friend. 
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BIOGRAPHICAL SKETCH OF PROFESSOR S.K.SRINIVASAN 



The following brief account of the life of Professor S. K. Srinivasan might 
help the reader to put him and his work in proper perspective. 

Professor Srinivasan was born in Kanchipurayn in Tamil Nadu, India on 
December 16, 1930. After completing his high school education in Cheyyar , he 
joined the Vivekananda college to study Mathematics and obtained his B. A. 
(Hons) degree of the University of Madras in 1953. Then he had to decide what 
to do next. This was a turning point for him and his inclination towards an 
academic career took shape then. The University of Madras had just started a 
new Department of Physics and attracted by the opportunity to work with 
Professor Alladi Ramakrishnan, he did his Master’s programme in Theoretical 
Physics there, receiving his M. Sc. degree (by research) in 1955. Subsequently 
he continued for the Ph. D. programme with Professor Alladi Ramakrishnan as his 
guide and submitted his thesis entitled "Theory of random integrals and their 
applications to physical problems" in 1957 and was awarded the Ph.D. degree of 
the University of Madras in 1958. Just as his later contributions, his thesis 
work had a fundamental character and was directed towards physical 
applications. 

After submitting his doctoral thesis, he went to the University of Sydney 
as a postdoctoral fellow, where he had the opportunity to work with Professor 
Messel and to familiarize himself with the first computer in Australia — 
SILLIAC. After a year he returned to the University of Madras, accepting a 
senior research fellowship of the National Institute of Sciences of India. 
Subsequently he joined the Indian Institute of Technology, Madras in 1959 at 
its very inception. Here he went on to become a professor in 1967 and a senior 
professor in 1974. We can only hint here at the magnitude of Professor 
Srinivasan’ s contribution to the Indian Institute of Technology, Madras in 
general and to the Department of Mathematics in particular. He spent a great 
deal of his time to improve and maintain the quality of the Department. This 
is true not only when he was the Head but also before and after. 

Teaching ability should be measured by one’s effectiveness as a teacher 
as well as one’s contribution to course design, curriculum development and 
continued improvement of teaching methods. Judged from these aspects 
Professor Srinivasan is a very good teacher. He was particularly creative with 
regard to the introduction of modern subjects. 
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At the Indian Institute of Technology, Madras, Professor Srinivasan has 
served as the Head of the Department of Mathematics, as the Dean of Research 
Programmes and as a Member of the Board of Governors. Each position demanded a 
different kind of management and administrative skill and he was quite 
successful in each of these positions. In addition he has served in a host of 
other important committees. 

Research competence should be measured on the basis of distinguished and 
continuing contribution to scholarship; contribution to scholarship is 
expressed mainly in the form of publications as well as evidence of ability to 
provide academic leadership. Judged by these norms Professor Srinivasan is an 
excellent researcher. In fact, the phenomenal growth of advanced research in 
the field of stochastic processes and their applications at the Indian 
Institute of Technology, Madras during the last three decades can be 
Justifiably attributed to the dynamic and purposeful stewardship of Professor 
Srinivasan. He has research interests in a large number of allied areas. To 
his credit he has 10 books and 170 research papers in established Journals. He 
has successfully guided 20 students in their Ph.D. programme. One of them got 
his degree in Electrical Engineering, one in Civil Engineering, two in Physics 
and the rest in Mathematics. 

No wonder several honours and distinctions came to Professor Srinivasan 
in the course of his brilliant and illustrious career as a scientist. Many of 
his publications have had wide recognition and have been quoted by several 
well-known authors. Professor Srinivasan has had visiting assignments at 
various Universities abroad: University of Sydney, State University of New 
York at Buffalo, University of Waterloo, National University of Singapore, 
Technical University of Munich and University of Texas at Austin. Professor 
Srinivasan is a life member of the Indian Mathematical Society, Operational 
Research Society of India (President, 1975-1976) and Bernoulli Society of 
Probability and Mathematical Statistics, an elected Fellow of the Indian 
Academy of Sciences and a founder Fellow of the Tamil Nadu Academy of 
Sciences. He was on the editorial board of "OPSEARCH" (1973-1977) and is a 
founder editor of the "Journal of Mathematical and Physical Sciences" now in 
its 24th year of publication. He is on the editorial board of "Solid Mechanics 
Archives" and "Annals of the Institute of Statistical Mathematics" and is a 
member of the Research Board of Advisers of the American Biographical 
Institute. Professor Srinivasan is listed in several directories of men of 
achievement and distinction. 
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His background in Physics combined with mathematical insight enables him 
to understand and relate to people interested in the applications of 
stochastic processes — probabi 1 ists, engineers, physicists, operations 
researchers etc. — and to bridge the gap among them. Professor Srinivasan 
feels completely at home in this diverse community. 

Professor Srinivasan is a original thinker and his work is always of the 
highest quality. His papers are marked by their clarity and originality. He is 
one of the early workers in the field of stochastic processes and to him are 
due some important results. A characteristic of Professor Srinivasan is the 
strong application bias, maintaining at the same time strong theoretical 
support and mathematical rigour for his work. His scientific curiosity and 
creativity seem to increase rather than decline over the years. In addition he 
was and continues to be tremendously helpful in advising and influencing many 
graduate students. He is selfless and extraordinarily generous in giving his 
time and concern to his students and colleagues. 

Professor Srinivasan is one of the most intellectual persons I have come 
across. I do not know anyone else who has as many interests as he and pursues 
them so deeply. He reads books like others read news papers. He is a great 
lover of music and likes gardening. He is a devoted walker. He is very 
pleasant and amiable. Any one who came into contact with him cannot but recall 
his acts of kindness and love. His warmth and affection reach out to any one 
who needs them. I am touched by the feelings expressed by a number of people. 
Discussions with them revealed a person who sometimes differed from the one I 
know. It is probably true that no one ever thoroughly knows another person in 
all his aspects and it made me feel good that there are other aspects to 
Professor Srinivasan* s personality that others had seen. He is truly a person 
worth knowing. 

It is very difficult to express how I admire him. The qualities of 
devotion, generosity, honesty and kindness are actualized in him, not from 
time to time but consistently in his daily life. He always sets the highest 
standards for himself, his collaborators and his students. He is an 
outstanding scholar and teacher. His career is a shining example of what a 
sincere and dedicated academician can achieve in this country in spite of the 
insurmountable obstacles in the path of progress. 



R. Subramanian . 




SCIENTIFIC CONTRIBUTIONS OF PROFESSOR S.K.SRINIVASAN 



All of Professor Srinivasan's work is in the area of Applied 
Probability. Even so, it is rather difficult to summarize his contributions, 
as they are scattered in different areas like the theory of stochastic 
processes and their applications to physics, operation research, biology and 
engineering. His work in the fifties is a consequence of his intimate contact 
with Professor Alladi Ramakrishnan and resulted in his doctoral dissertation 
on the theory and applications of stochastic processes in physics and allied 
fields which in turn laid strong foundations for his work in cascade theory, 
stochastic integrals and differential equations dealing with noise phenomena. 
His familiarity with cascade shower theory naturally provided the motivation 
and desire for his investigations in high energy physics, pion physics and 
scattering theory. While he had adequate reasons for pursuing theoretical 
physics, with the appointment to a full professorship in applied mathematics 
at the Indian Institute of Technology, Madras in 1967, he had shifted his 
interests to other typical applied probability areas like operations research 
and biology, besides noise and fluctuating phenomena where there was abundant 
scope for the study of specific non-Markov processes. 

Theory of Point Processes: 

Professor Srinivasan* s contributions in the area of point processes are 
manifold - extension of the theory of product densities to the 
multidimensional and non-Euclidean spaces on the one hand and irregular point 
processes on the other. His doctoral work on cascade theory was developed 
further in the early sixties in collaboration with Professor Iyer, 
Dr.Koteswara Rao and Professor Vasudevan and this led to the evolutionary 
sequent correlation functions as a means for the study of the general theory 
of point processes. His early research work on twins and multiplets in the age 
dependent population growth culminated in 1961 in the formal theory of 
multiple point processes and multiple product densities. Apart from this he 
had introduced several innovations for the study of special point processes 
and in particular he characterized the space charge limited shot noise 
processes (now studied in the literature as inhibited point processes) and the 
resulting cumulative response processes. A significant outcome of his 
investigations is the emergence of the new formalism for the study of the 
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kinetic theory, fluctuation phenomena and onset of turbulence apart from 
original ideas on second order point processes. The formulation of stochastic 
integrals and differential equations through the sample path analysis engaged 
his attention in the fifties and the encouragement of Professor Richard 
Bellman and his teacher Professor Alladi Ramakrishnan enabled him to take up a 
systematic study of the processes arising in physics and engineering 
(Srinivasan and Vasudevan, 1971). His early work in this relates to the 
characterization of processes arising from differential equations where the 
forcing terms are stochastic in nature (Srinivasan and Mathews , 1956, 
Srinivasan, 1960). With the advent of McShane-Ito calculus he was naturally 
led to many innovations culminating in an appropriate definition of stochastic 
differentials of jump processes, more specifically of point processes 
(Srinivasan, 1977, 1982). 

In the late sixties he started investigating interacting point processes 
with Dr . Ra jamannar and Dr . Rangan and made i nno vat i ons i n mode lling firing 
sequence of single neurons. The interacting point processes or for that matter 
even renewal processes are a bit tricky to handle; his main contribution lies 
in the discovery of the appropriate regenerative structure which in turn 
enables one to obtain the characteristics of the process. These are codified 
in his work on single neurons. 

Cosmic Rays and Elementary Particle Physics: 

The cascade theory of cosmic ray showers as formulated by Bhabhaand 
Heitler deals with the distribution of particles over energy at a particular 
thickness. Motivated by the experimental findings in the mid-fifties, 
Professor Srinivasan in collaboration with Professor Ramakrishnan came up with 
a modified approach and expressed the expected value of the size of th > shower 
in terms of production product densities defined over a two dimensional 
continuum. The estimation of the fluctuation of the size, which was normally 
considered to be formidably difficult to compute, was carried out in 
collaboration with Professor Iyer and Dr.Koteswara Rao (Srinivasan, 1961, 
Srinivasan and Iyer and Koteswara rao, 1964) through a set of compact 
equations obtained by the use of the invariant imbedding technique. This work 
on high energy cosmic rays has stimulated him further to take up 
investigations in high energy physics where he made significant contributions 
in the area of multiple production of particles. Some of his recent 
contributions in the area of multiplicity distribution and quark gluon cascade 
are further innovations. 
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Operations Research: 

In the sixties he was attracted to the area of operations research due to 
the rich non-Markov structures present in the stochastic models of queues, 
inventories and reliability. Partly independently and partly in collaboration 
with Professor Subramanian and Professor Kalpakam, he developed a unified 
approach by discovering regenerative structures in the stochastic processes 
governing these models. We wish to make a special mention of the complete 
solution for the most general type of (s,S) inventory model (Srinivasan, 1979) 
and the analysis of redundant systems (Srinivasan and Subramanian, 1980). All 
along Professor Srinivasan had the conviction that the study of point 
processes would lead to the characterization of non-Markov processes. Here he 
was able to satisfy himself and also convince fellow workers in ample measure. 

Optics: 

In the area of optics, he was drawn towards the statistical detection 
theory. His work in the sixties relates to the correlational structure of the 
point process of photo electron emissions of partially coherent light 
(Srinivasan and Vasudevan, 1971) . He then undertook detailed investigations 
leading to the determination of photo electron statistics of light beam with 
arbitrary spectral profile (Srinivasan and Sukavanam, 1971, Sudarshan et ai . , 
1973, Srinivasan, 1974). During the past seven years he is working on the 
population theoretic approach to the problem of light and its detection; as a 
result of his recent findings we now have a complete population point 
processes theory of light and its detection. An important outcome of this 
investigation is the identification of a class of doubly stochastic Poisson 
processes with the point process generated by the emigrations of a population 
process (Srinivasan, 1988). His latest work relating to intermittency in 
particle production and his projected work on quantum optics will introduce 
new vistas in the field. 



R . Subramanian . 
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1. Introduction 

Cascade theory of cosmic ray showers has provided very good motivation 
for the development of the general theory of stochastic processes. Cosmic ray 
showers have been observed around 1930 when stochastic processes have not 
been studied in depth. The experimental data observed around that period 
needed a comprehensive stochastic model to comprise in itself the twin aspect 
of branching with respect to one parameter (energy) and evolution with respect 
to a different parameter (thickness of absorber). Unfortunately the study of 
branching processes was taken up by probabi lists only much later; the 
celebrated work of Kolmogorov and Dmitriev leading to the limiting behaviour 
of population size emerged around 1947. The early work on point processes 
dates around 1943 when Palm published his work on streams of calls in a 
telephone exchange. Although population theory can be traced to the time of 
Galton and Watson, stochastic problems of population growth received attention 
only during forties. Thus when the experimental data on cosmic ray showers 
confronted the people around thirties, there was a search for a possible 
interpretation in terms of a statistical theory of population of points with 
specific labels. In 1937 Bhabha and He i tier and Carlson and Oppenheimer 
independently came out with their formulation in terms of what are now known 
as inclusive probability measures and thus laid a strong foundation for the 
theory of population point processes. Bhabha soon realized that his 
formulation of cascade showers was only a first step and was on the quest for 
further tools for the determination of the fluctuation in the size of the 
shower about its average value. His efforts culminated in a paper (jointly 
authored with Ramakrishnan) which provided a frame work for the calculation of 
the size fluctuation of showers. The work, from the probabi 1 ist’ s point of 
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view, is highly significant in as much as it provided an explicit example of a 
bivariate population point process which is evolutionary rather than 
stationary in character. Both Bhabha and Ramakrishnan realized the importance 
and relevance of the work to the general theory of stochastic processes and 
soon independently published their theory of point distributions. By the time 
these results were published, there had been significant developments in the 
theory of population growth partly due to axiomat izat ion by Kolmogorov and 
partly due to investigations in Galton-Watson processes. David Kendall had 
just then introduced his cumulant functions to deal with age structure in 
stochastic population growth. Quite unaware of Palm’s contributions, Bellman 
and Harris developed the method of regeneration point by dealing with 
age-dependent population growth and were ultimately led to versions of limit 
theorems on population size of a nature more general than those attempted by 
Kolmogorov and Dmitriev (1947). After the publication of the results by 
Bhabha and Ramakrishnan, Bartlett and Kendall (1951) reformulated the problems 
in terms of differential equations satisfied by the characteristic functional 
of the population which has meanwhile been introduced by LeCam and Bochner. 
By an elegant application of the regeneration point technique , they were able 
to derive the differential equation satisfied by the characteristic functional 
which in turn yielded the Bhabha- Ramakrishnan equations in the context of 
shower theory on the one hand and the Bellman Harris age-dependent population 
growth equations on the other. 

The story of cascade theory of showers is a long and interesting one to 
narrate; it cuts across many areas of the theory of stochastic processes and 
is highly suggestive of many other applications in other disciplines as well. 
The theory of population point processes which was the main offshoot remained 
to be an area of fruitful research until quite recently. Hence I feel it is 
worthwhile to present, in proper prospective, a brief account of the manner in 
which cascade theory of cosmic ray showers evolved and led to a\ rich theory 
of population point process. 

2. Electron Photon Showers 

To appreciate the impact of cascade theory on the development of the 
theory of stochastic processes, it is necessary to recall some of the salient 
features of the (experimental) cosmic ray research around 1912 when the extra 
terrestrial origin of radiation in the atmosphere was inferred from the baloon 
experiments of Victor Hess. It was soon demonstrated that cosmic rays at 
ground level can initiate cascades and are stopped in lead plates. With the 
advent of quantum mechanics, the cascade showers of particles generated by 
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electrons or photons were sought to be explained by means of the two 
fundamental processes: Bremmst rah lung (radiation) by electrons and pair 

production by photons. Bhabha and Heitler (1937) analyzed the showers 
generated by a single primary by dealing with the expected value of the number 
of electrons (each with energy > E) that are detected at a definite depth of 
atmosphere (measured in an appropriate unit called cascade unit). By varying 
the parameter t representing the depth, they wrote down differential equations 
representing the differential mean (expected value) number of electrons and 
photons with energies in (E, E + dE). Carlson and Oppenheimer (1937) 
simultaneously and independently dealt with the problem using an identical 
approach. The main conclusion was that the mean size of the observed showers 
could be explained in a statistical frame work on the basis of the two 
fundamental processes: Bremmst rah lung and pair production. The problem also 
attracted the attention of the Russian Physicist Landau who in collaboration 
with Rumer came out with his own analysis providing some insight into energy 
degradation of individual electrons. In current day terminology, Landau’s 
treatment is essentially a sample path analysis of continuous parameter Markov 
process. To appreciate this we note a typical shower can be represented as in 
Figure 1 .Solid lines represent electrons while wavy lines represent photons. 
The problem in the first instance was treated as one dimensional (with 
reference to depth of atmosphere) the lateral spread being neglected. If we 
persue the primary electron, its life history is easily narrated; it goes on 
emitting photons (due to interaction with matter) in a probabilistic sense 
until its energy is reduced considerably that it cannot emit photons of 
substantial energy. It should be appreciated that Landau and Rumer recognized 
the life history of an individual electron as a Markov process with respect 
to t, the state space being the one dimensional continuum representing energy. 
A rigorous mathematical treatment was later provided by Harris (1957). 

3. Point Distributions and Product Densities 

To get a proper perspective, we introduce the notation N^E.t) 

representing the number of particles of type i ( electrons correspond to i=l 

* 

while photons to i=2) with energy < E at a definite thickness t . The 

function N^Ejt) is a random variable for a fixed t. At t = 0, we have a 

single electron of energy E . The quantum mechanical processes imply that the 

o 

electron with energy E has a probability R (E,E*)dE’ per unit thickness (of 
material traversed) of radiating a photon itself dropping to an energy in the 
interval (E’,E’+dE’). Likewise a photon (if sufficiently energetic)of energy 
E has a probability R 2 (E,E’)dE’ per unit thickness of producing an 
electron-positron pair one of which has an energy in the interval (E’ ,E’+dE’). 




XXXI 



Thus it is clear that the physical processes imply that an incident electron 
or photon gives rise to a shower as in Figure 1. The main problem was how to 
arrive at a frame work within which the statistical characteristics of the 
shower of particles can be described. For each t, we have particles 
(electrons/photons) distributed over a continuum (energy). Bhabha and Heitler 
dealt with the expected value of N. (E, t) and wrote down the conservation 
equation w. r. t . t. 

The year* 1937 appears to be a crucial one; in kinetic theory Yvon, a 
French Physicist, (see also de Boer) had to deal with a similar problem of 
distribution of a discrete population over a continuous parameter. Here the 
parameter represents the phase space and the population collection of 
molecules. Yvon not only dealt with the expected value of sum functions 
associated with the position momentum of the molecules but also its mean 
square value through an appropriate second order correlation function. In 
modern terminology of point processes, Yvon’s formula corresponds to the mean 
square value of cumulative response of an arbitrary weighted point process. 
Apparently neither Bhabha nor any of the other proponents of shower theory was 
aware of Yvon’s contribution. Furry and Arley (1943) studied the problem and 
they were able to propose simple multiplicative models starting from Poisson 
and population birth processes. Scott and Uhlenbeck (1942) circumvented the 
difficult problem (of describing the distribution of a discrete population 
over a continuum) by dividing the energy interval into a finite number of 
portions and dealing with the population in each of the portions; however they 
made some approximations which eliminated the correlation that otherwise would 
have persisted. Bhabha proceeded to refine his calculations by taking into 
account loss in energy by electrons due to ionization of the atoms in the 
media; in collaboration with Chakrabarti (1943, 48), he obtained the explicit 
Mel 1 in transform solution and numerical estimates for the average number 
(first moment) of electrons. Since there was no awareness of Yvon’s 
contribution on the part of shower theorists (and that too despite Bogliubov’s 
reformulation of kinetic theory in 1946), the problem of estimation of the 
fluctuation of the size of the showers about its mean value remained open. 
Thus Bhabha continued his search for a new mathematical technique that would 
yield a formula for the mean square value of the number of electrons and 
ultimately succeeded in 1950 by building up a theory of point distribution. 
He essentially viewed the inclusive distribution which formed the basis for 
the differential mean numbers of electrons as the one arising from appropriate 
summation/ integration over the points of the finite dimensional distribution. 
Bhabha dealt with a continuous parametric assembly of particles; he considered 
a possible eigenstate in which there are m systems (particles) present with 
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the parametric values in the intervals (x ,x +dx ) (x ,x +dx ) ... 

111222 

(x ,x +dx ), (x <x <x <x ) and introduced the probability measure 

m m m 12 3 m 

\fi (x ,x . ...,x )dx dx ....dx with the help of which he was able to obtain the 

ml2 m 12 m 

expected value of any arbitrary function of the parameter. By introducing 
appropriate set functions he was able to obtain moment formulae. As an 
application of the formula, he in col loborat ion with Ramakrishnan, dealt with 
the problem of electron-photon cascade and wrote down the extension of the 
diffusion equations for the second order functions expressing the correlation 
in energy (for any fixed t). Since there are two types of particles, this 
gave rise to four functions leading to a linear simultaneous system of integro 
differential equations. By the combined use of matrix approach and Mel 1 in 
transform, they were able to obtain an explicit solution and finally the 
second moment was expressed as double Laplace inversion integral. Although 
the expression was unwieldy, it was amenable to numerical inversion. In 1954, 
Ramakrishnan and Mathews prepared tables of variances of the size of the 
showers for various primary secondary energy ratios. 

The probability measure induced by the density \fj can give rise to 

m 

inclusive probability measures. As this is an important point, it is worth 
discussing the same in some detail. If we fix the parameter x^and integrate 
over the rest of the parametric values x ... x taking into account the 

2 m 

indistinguishability of the systems except by their parametric labels, we 
obtain 



f>) = T — 1 1 ^ [ f. . . f 0 (x, x .... x ) dx . ..dx 
11 (.m-lj!JJ J m 1 2 m 2 in 

where the integration is over the entire range of the parametric values. Now 

the quantity f^Cx^dx^ has a simple interpretation in as much it denotes the 

probability of finding m systems with a typical member having parametric 

value in (x , x + dx ) . If at this stage we sum over all values of m 
111 

running from 1 to 00 and denote the resulting function by f (x ), then 

f i (x^)dx i denotes the inclusive probability of finding a system in 

(x i ,x i +dx i ). It is precisely this function for which differential equations 

were provided by Bhabha and Heitler and Carlson and Oppenhiemer On further 

integration of x^ over any specific range, we obtain the expected value of the 

number of such systems in that range. An extension of this argument naturally 

leads to correlation functions of the type f (x , x , ... x ) where 

n 1 2 n 

f (x , x , ...x ) dx dx ... dx denotes the probability of finding 

simultaneously systems in parametric ranges (x. , x^ + dx^ ) (i = 1, 2, ... n) 

irrespective of the number found elsewhere. Bhabha attempted to express the 

moment formulae in terms of these functions. Ramakrishnan, independently and 

around the same time, introduced the functions f directly through the random 

n 
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variable N(x) representing the number of particles in the interval [0, x] and 
obtained an explicit formula for the moments of the number of particles with 
parametric values in an arbitrary interval [a, b]. He called these functions 
as product densities; in fact such functions were introduced earlier by Rice 
(1945) in the context of zero crossings of a Gaussian process and were known 
under the name inclusive probability densities. Apparently such functions 
have been in vogue in kinetic theory of fluids (see Bogoliubov (1946)). 
However the merit of the work of Ramakrishnan lies in the moment formula and 
the formulation of evolution equations for the product densities. 

The probability measures induced by the functions 0 were also studied by 

m 

Janossy quite independently in the same year 1950. Janossy in fact studied 
the problem of showers and being aware of the inherent difficulty due to lack 
of general methods of handling such point distributions introduced in 1949 
his famous G- equation by dealing with the mass function n (n,E,E ,t) 

i o 

representing the probability that there are n electrons at depth t each with 

energy > E in a shower generated by a primary of type i and energy E . 

o 

Janossy used the fact that once a secondary particle is produced by the 
primary, the shower generated by the secondary is statistically independent of 
the residual shower. This in turn enabled him to obtain an equation, although 
integral in nature, for the ?r-f unctions: 

t -p (E )x 

e Y n (n ,E,E’ ,t-x) 

La 11 



-p .(e ) t 

n (n ,E,E -E’,t-x) R (E\ E ) dE dx + S e 

3“i 2 o i o il 

where 

E 

r ° 

p (E ) = R(E' I E )dE' 

i o J i 1 o 

If we now define 



G.Cu.E.E ,t) = V u n it (n,E,E ,t) 

i o i o 



we obtain 

n t —p ( E ) x 

e ° G^u, E, E’ , t-x) 

p (E )t 

G (u,E,E -E\t-x)R (E’,E }dE’ dx + u S e 

3~i o i o il 
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Janossy’ s idea stepped up the pace of research and one could calculate at 
least in principle the various moments of the electron distribution. In fact 
the equation itself is a generalization of the Bellman-Harris equation 
proposed in the context of age dependent population growth. In the same year 
Janossy came up with the next contribution wherein he dealt with the iff 

m 

functions and indicated how diffusion equations for iff can be written down; 

m 

he also indicated a method of obtaining the moments of the number of particles 
produced in a shower upto a depth t. Messel and his collaborators have 
always preferred to use the Janossy measure as the basis in as much as it is 
easy to justify the use of forward differential equation as a consequence of 
the Markov property being imposed on Janossy measure functions. The impact of 
the work of Janossy can be judged from the fact that more than fifty papers on 
the subject particularly by Messel and his collaborators appeared in the next 
few years. 

It is rather difficult to accord priorities In these contributions. Point 
distributions are perhaps known to theoretical physicists well before 1950; in 
particular the inclusive probabilities have been used in kinetic theory as 
early as 1937. Bhabha used the iff functions to arrive at the inclusive 
distribution; Ramakrishnan directly dealt with the inclusive distribution. 
Although diffusion equations were arrived at using physical arguments of 
number conservation in the first instance, a formal proof based on the Markov 
property enjoyed by the Janossy measures was provided a year later 
(Ramakrishnan 1951). Janossy introduced his measure quite independently and 
later on used Markov property to arrive at the equations for inclusive density 
functions. This contribution has been missed in view of his more important 
contribution leading to the celebrated G- equation. 

4. Further Developments in Shower Theory 

The second moment of the number of electrons in a typical electron photon 
shower as calculated by Bhabha - Ramakrishnan or Janossy and Messel was more 
or less in broad agreement with the experimental data and it was concluded 
that except for same highly energetic primary where there was a possibility of 
direct pair production of electron positron by electrons, the data more or 
less gave conclusive evidence of the showers being a consequence of the two 
fundamental processes of quantum mechanics. However around 1954 there were 
anamolous showers and a more extensive investigation particularly by Fay at 
Gottingen was taken up wherein showers were observed in emulsions; it was 
pointed out that it would be more convenient to count the particles with 
energy labels at the points of production rather than at any later depth. To 




XXXV 



facilitate comparison with experimental data, Ramakrishnan and Srinivasan came 
up with a new approach wherein the deal with electrons produced over a certain 
thickness say an interval [0,t] with energy of each particle (or pairs that 
are observed) is greater than E at the point of production. Thus in the 
product density function f (x , x , ... x ) the parametric space was 

n 1 2 n 

identified to be the Cartesian product of E (energy) and t (depth). In this 
case, the main quantity of interest is the expected number of electrons 
produced in a small stack of emulsions. Tables of mean numbers were prepared 
for the various energy ranges using the then powerful computer SILLIAC and the 
analysis done by Srinivasan in collaboration with Messel and others (1958) 
showed that the anamalous showers could be explained within the framework of 
the basic quantum mechanical processes. The investigations gave a new lease 
of life for the cascade theory for two reasons. First the problem naturally 
led to a two dimensional point process. Second the calculation of higher 
moments had to be taken up so that the main conclusion that there was no 
anamaly could be confirmed beyond doubt. This time there was an added 
advantage. Janossy’s G- equation could be adapted to suit the special 
requirements of the problem. Srinivasan (1961) came up with an idea that the 
unwieldy Bhabha - Ramakrishnan equations could be given up and in their place 
a new set of compact equations for the moments smaller in dimension arrived at 
by a judicious use of Janossy’s regeneration point method. These equations 
were then solved by recourse to Mel 1 in transform technique and by standard 
complex variable methods, inversion of the expressions were executed elegantly 
to yield simple formulae for the moments upto 3rd order. (Srinivasan, Iyer 
and Koteswara Rao (1964)). 

There was also a parallel development in the theory of showers 
particularly by Dyson and Me Voy who were motivated by some of the existing 
results in parity non-conservation and were led to the investigation of high 
energy showers initiated by electrons with a specific spin. There arose an 
interesting question whether the polarization properties were transferred to 
the secondary particles in the showers. Dyson and Me Voy proceeded on the 
same lines as Carlson and Oppenhiemer and answered in the affirmative by 
dealing with the expected value of the number of electrons with polarization 
and energy labels. Koteswara Rao in his doctoral dissertation (1967) carried 
out the analysis and provided compact formulae for the various higher moments 
arising from such population point processes. 

5. Population Point Processes and Further Outlook 



The above developments in physics have naturally provided ample motivation 
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for the construction of a formal mathematical theory of population point 
processes. The theory of product densities was generalized to include within 
its frame work multiple points or points which are themselves populations. 
For the cascade theory itself other questions like the non-finite nature of 
the quantum mechanical cross-sections were also resolved; the monograph of 
Harris (1963) contains an excellent summary of the salient features of 
electron- phot on cascades as well as point distributions in general. With the 
publication of a monograph on analytic methods (Srinivasan 1969) leading to 
compact moment structures, it was generally believed that research in the area 
of cascade theory had reached a satisfactory state of culmination and that the 
techniques presented therein would be used in contexts other than cascade 
theory. However this is not to be; while there has been a wide appreciation 
of population point processes in areas like statistical physics, 
neuro-physiology, light amplification and detection (see for example Uhlenbeck 
(1971), Murthy (1974), Holden (1976), Sampath and Srinivasan (1977), Shephered 
(1981, 1987), Teich et al (1984) Srinivasan (1988)), there is also a marked 
revival of activity in the theory of cascade processes with special reference 
to characteristics like energy and rapidity in jets observed in collider 
experiments. Some of these activities are spearheaded by Carruthers and Shih 
(1987), Giovannini (1972, 79) and Van Hove (1986). The Bellman Harris 
regeneration equation and the improvized imbedding equation again due to 
Bellman, Kalaba and Wing (1960) have reappeared in the form of 
Altreilli-Parisi equation (See Hwa ( 1988)). Right through seventies, Feynman 
in collaboration with Field (1978) (See also Fukuda and Iso (1977)) studied 
the jets and being unaware of some of the powerful analytic tools advocated 
extensive use of Marte Carlo simulation techniques. In fact the Feymman-Field 
formulation of e + e collisions resulting in big jets is essentially a 
reformulation of Janossy’s nucleon cascade model in terms of the production 
product densities introduced earlier by Ramakrishnan and Srinivasan (1956) for 
the study of QED cascades in emulsions. In many workshops, there are tutorial 
sessions on branching models mainly to understand the particle multiplicity 
distributions in jets produced in colliders. Attempts are also being made to 
examine the applicability of some of the distributions used in light beam 
detection (Fowler et al 1988, Voordas and Weiner 1988); in this connection it 
is worth noting that light amplification and detection process is by itself an 
interesting population point process. 

There are also important attempts to use the ideas that are currently in 
vogue in the modern theory of turbulence particularly like self-similarity, 
intermittency and fractal nature. (See for example Bialas et al (1988), Chiu 
and Hwa (1990)). Since the new particles and their interactions by themselves 
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constitute a complex system, currently some investigations are on to verify 

whether the electron-photon cascade showers do enjoy some of these properties. 

These investigations are really opening up new vistas in our understanding. 

Let me illustrate by a simple problem which I have examined in the last few 

months in collaboration with Chiu of UT at Austin. In the cascade theory the 

main object of interest is the distribution of the number of electrons 

produced with energies > E . To examine the existence of fractal structure in 

o 

the distribution or self-similarity in energy ranges, it is natural to look at 
the electrons with energy in a window (€,€+5). There are two methods of 
analysis available : We can introduce an innovation in the Janossy G-equation 

by introducing another label to characterize the configuration in which the 
primary is found, since the primary can be to the left or right of the window 
or in the window itself. The G- equation can be written down and it in turn 
yields recursive integral equations for the moments. Right now we have found 
that although the equations have attractive features, quick inferences are not 
possible. Hence we examined the possibility of dealing with the product 
densities of the cascade process; it has turned out to our surprise that it is 
possible to obtain explicit analytic expressions for the second order product 
densities which in turn will enable us to obtain the moment structure of the 
number of particles in the window at least to the second order. The 
investigations are still continuing and the findings will enable us to 
conclude whether the branching nature will give rise to intermittency and 
fractal structure in multiplicity distribution. 

In summary you will agree with me if I say that we have had an exciting 
period of over fifty years of research in cascade theory of showers that had 
dominated our activities of research in the area of stochastic processes and 
more particularly in population point processes. 
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Figure 1: A typical shower of electrons and photons 

produced by a primary electron. 
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ABSTRACT 

The power spectral density of a stationary square wave process with sign 
changes at the transitions of a finite state Markov renewal process is 
derived. Particularly explicit formulas are obtained for the Markovian 
Arrival Process, a generalization of the Poisson process with a natural matrix 
formalism which commonly leads to useful explicit results. An application to 
a procedure for the quantification of burst i ness is discussed. 

1. INTRODUCTION 

This paper deals with a descriptor of the random point process generated 
by the transitions of a finite-state Markov renewal process. We begin by 
defining some essential notation. We consider am irreducible, positive 
recurrent, m-state Markov renewal process with transition probability matrix 
H(.). The matrix Laplace-Stielt jes transform of H(.) is denoted by h(s) and 
the finite positive vector of row sum means -h' (0+)e by £. The column vector 
e has all its components equal to one. The stochastic matrix H = H(oo) has the 
invariant probability vector n. The matrix A(£) is a diagonal matrix with the 
components of £ as its diagonal elements. The inner product E = t r|S is called 
the fundamental mean of the Markov renewal process. It plays an important 
role in many asymptotic results. The quantity X = E *, the fundamental 
rate , is the rate at which transitions occur in the stationary version of the 
Markov renewal process. 

It is well-known (see Pyke [13] or Cinlar [4]) that the stationary 
version of the Markov renewal process is obtained by choosing the initial 
state according to the probability vector E ~ 1 n A(0) and by making the first 
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transition according to the semi-Markov matrix H*(x), defined by 

H (x) = A [H - H(u)] du, for x £ 0 . 

In this paper we shall only deal with the stationary Markov renewal process so 

constructed. In order to avoid uninteresting issues, we shall also assume 
that, with probability one, there are only single transitions. That 
assumption is entirely "technical" and is satisfied in most applications. 

For a stationary point process with single transitions, we may define an 
associated square wave process. We consider the counting process which is the 
random counting measure N induced by the stationary process on the Borel 

setsof R. In what follows, it is sufficient to consider the random variables 

N{ [0, t ] } = N(t) = inf {n : S s t>, for t * 0 . 

n 

The random square wave corresponding to the stationary point process is 
the stochastic process defined by 

Y^O) = 1 or 0, each with probability 1/2, 

and 

Yl (t) = max <0,(-l)V o)+N<[0 ’ tl>+1 > , for t real . 

The process {Y^.)} is clearly stationary and E [Y^t)] =1/2 for all t. The 
centered square wave is simply the process Y(t) = Y^t) - 1/2 . 

The square wave power spectral density (SQSD) of the stationary point 
process is defined by 



00 

S(f) = [ e‘ lwt R(t)dt (1) 

-00 

where u> = 2nf and R(t) is the autocovariance function 

R(t ) = E [Y(0) Y( t )] 

of the process Y(.). The function R(.) is clearly even, and for t £0, we have 

R(t) = E [Y^Oiyt)] - i i P{N(t) is even > - j (2) 

The choice of the variable f, representing the "frequency" in cycles rather 
than in radians, is more prevalent in the engineering literature. Whether v> 
or f is used is entirely a matter of convention. 
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In terms of the probability generating function 

P*(z,t) = E {z N(t) >, for t £ 0 (3) 

we may rewrite (2) as 

R(t) = i P*(-l,t) (4) 

The preceding considerations are well-known, at least for the Poisson 
process [5). For a Poisson process of rate A, one easily obtains that 

R(t) = ^ exp {— 2A 1 1 1 > , for all real t (5) 

so that the SQSD is given by the Cauchy "density" 

S(f ) = J [X 2 +It 2 f 2 ]’ 1 (6) 

For other instances of uses or discussions of the SQSD, see Castro, Kemperman 
and Trabka [3], Lamond [6,7], Neuts and Sitaraman [14] and Neuts [12]. In 
this paper, we shall give a derivation of a matrix-analytic expression for the 
SQSD of the point process generated by the transition epochs of a finite state 
Markov renewal process. We shall discuss an important particular case for 
which that density is highly tractable and conclude by describing a potential 
practical use of the SQSD . 

2. THE DISTRIBUTION OF THE COUNTS 

The derivation of the probability distribution of the counting variable 
N(t) proceeds in classical fashion as in Pyke [13]. The following expressions 
are given for completeness. For the stationary Markov renewal process, we 
define 

00 

P(n; t) = P{N(t) = n| N(0) =0>, and P*(n;s) = j e" st P(n;t)dt, for n £ 0 . 

O 

By direct probability arguments, we obtain that 

t 

P(0; t) = E'Sr J [H - H(u)]du, and P*(0;s) = s^-s'V 1 * [I-h(s)]e (7) 

O 

and for n £ 1 , 

r l r 1 ” 11 t i ^ 

P(n;t) = E'V I [H-H(u)]du I dH (n '(v)[H-H(t-u-v)]e (8) 



so that 
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P*(n;s) = s"" 2 E"" 1 7i [I -h(s)] 2 h n 1 (s)e (9) 

Routine calculations lead to the generating function 

P (z;s) = £ P (njs)z 11 : 

n- o 

P* ( z ; s ) = s'" 1 - (1 -z)s” 2 E“ 1 tt [I - h(s)] [I - zh(s)] 1 e (10) 

The Markovian Arrival Process : A particularly tractable class of Markov 

renewal processes are those whose transition probability matrix is given by 

r x 

H(x) = exp(Cu)du D, for x £ 0 (11) 

o 

where C and D are square matrices of order m whose sum Q = C + D is an 
irreducible infinitesimal generator. The matrix C has negative diagonal 
elements, nonnegative off-diagonal elements and its inverse exists. The 
matrix D is nonnegative. 

Such a Markov renewal process (and the extension allowing for group 
arrivals) is known as a Markovian arrival process , and because of its 
appealing matrix-analytic properties, has already been extensively, used in 
queueing theory. We refer to Asmussen and Ramaswami [1], Blondia [2], 

Lucantoni [8], Lucantoni, Meier-Hellstern and Neucs [9], Neuts [11], for 
detailed discussions. Particular Markovian arrival processes are the PH- 
renewal process and the Markov-modulated Poisson process [11] and 

(independent) superpositions of these . 

If we denote by © the stationary probability vector of the generator Q, 
then it is readily verified that X* = © D e, and n = (A ) 1 © D. The transform 
matrix h(s) is given by 

h(s) = CsI-C) -1 D (12) 

By substitution in (10) or by a direct argument as in Neuts [11], we verify 
that 

P*(z;t) = e(sI-C-zD) -1 e (13) 

For the sake of caution, we note that for some important cases, the stochastic 
matrix H = (-C) _1 D, is reducible (because some columns of D vanish, ) but that 
rarely affects the validity of the calculations based on the useful matrix 
formalism for the Markovian arrival process. It is more convenient to 
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conserve the matrix formalism than to work with the reduced transition 
probability matrices. 

3. THE SQUARE WAVE SPECTRAL DENSITY 

Theorem 1 : For a non- lattice, finite-state Markov renewal process, the 
SQSD is given by 

S(f) = u'Vi [I + h( it*>) ] _1 [I - hCica)h(-ico)] [I + hC-iw)]' 1 ® (14) 

for f ^ 0, and where cj = 2nf . 

Proof : By virtue o* (4), the Laplace transform R (s) of R(t) is given by 

R*(s) = i P*(-l,s) = L - K I [I-h(s)l [I+h(s)] _1 e (15) 

4 2s 2 E 

Now observing that the Fourier integral in (1) may be written as 

R (io>) + R (-iu>), the expression for S(f) is obtained by routine matrix 

manipulations. 

Remarks : By rather belabored but straightforward calculations, we may 
express S(f) in (14) in terms of the real and imaginary parts of the 
characteristic matrix h(iw) (whose complex conjugate is the matrix h(-io>)). 
Since a complex matrix and its conjugate do not commute in general, the "real" 
expression for S(f) is more involved than that given in (14). 

The requirement that H(.) is not a semi-Markov matrix of lattice type is 

essential. For that case, the square wave spectrum is discrete and a separate 
analysis is required. When H(.) is non-lattice, the inverses in (14) exist by 
virtue of a corollary to the Perron-Frobenius theorem. 

In what follows, we shall obtain an expression for the value of S(0+) by 

series expansions of the matrices involved in formula (14). In that analysis, 
the case where the stochastic matrix H has -1 as an eigenvalue, or 

equivalently where the embedded Markov chain of the Markov renewal process is 
periodic with an even period, requires special treatment. The calculations 
for that case are, in general, quite belabored but follow the same lines as in 
Lamond [ 7 ] . 

Henceforth, we assume that the transition probability matrix H( . ) is 

non-lattice and that all sojourn time distributions have finite second 
moments. We denote the first and second order moment matrices respectively by 

H = -h' (0+) and H =h"(0+) 

1 2 
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Theorem 2 : Provided that the matrix I+H is nonsingular and under the 
stated assumptions on the moments of H(.) f the value SCO) of the SQSD at zero 
is given by 

* 

S(0) = {7tH 2 ( I + H) _1 e - 2 tcH i ( I + H) _1 H (I+H) -1 e> (16) 

Proof : From the expansion 

2 

h(iw) = H - iwH - %- H + o(w 2 ) 

12 2 

is follows that 

2 

7r [I - h( io>) ] = iwrcH + ~ ttH 2 + o(w 2 ) 

Writing 

2 

[I + h(iw)] _1 e = v(0) + iwvCl) v(2) + oCw 2 ) 

we obtain, provided that I+H is nonsingular, the following expressions for the 
coefficient vectors v(0), v(l) and v(2) : 

v(0) = (I + H)‘ l e, and v(l) = (I + H) -l H (I + HJ^e 

v(2) = 2(1 + H) _1 H (I + H) _1 H (I + H)' 1 e - (I +H) _1 H (I +H) _1 e 
11 2 

Using these formulas, we compute the expansion of 

w [I - h(iw)] [I + hCiw)]” 1 © 

and add its complex conjugate to it. The terms in u> and w cancel. Upon 
substitution into expression for R* (io>) + R* (-iw),the stated result follows. 

For the Markovian Arrival Process , most of the preceding expressions 
greatly simplify. By virtue of formula (4) , 

R*( s ) = i e (si - C + D) _1 e 

and since e(C - D) = e(Q - 2D) =-2© D, we obtain by direct matrix calculations 
that 

S(f) = e D[4ir z f 2 I + (C - D) 2 ] -1 e (17) 

AN APPLICATION ; THE CLUSTERING OF A MARKOV RENEWAL PROCESS 

The qualitative features of the point process generated by the 
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transitions of a finite Markov renewal process can be very involved. We have 
been exploring analytic tools to quantify, for example, the time scales at 
which various degrees of clustering of the transition epochs manifest 
themselves. The SQSD is one such tool and the following discussion deals with 
its potential use in the qualitative analysis of such point processes. In 
that discussion, it should be borne in mind that a prevalence of "long gaps" 
in the point process will result in relatively higher spectral density values 
in the "low frequencies", whereas "dense clusters" tend to increase the 
spectral density at the higher frequencies. 

We can consider various forms of thinning which progressively "remove" 
clusters of points from the point process. By studying the SQSD for various 
values of a "thinning parameter" we can expect to gain insight in the levels 
of clustering in the original process. Specifically, a simple form of 
thinning consists of imagining a Poisson process of rate a , independent of 
the given point process. We agree that a point of the given process is 
"preserved" if and only if there is at least one Poisson event between it and 
its predecessor in the given process. 

The thinning corresponds to considering the arrival epochs of customers 
who receive service in a queue of capacity one, with a single exponential 
server and having the given Markov renewal process as its input process. At 
the expense of additional notation, we may readily replace the Poisson process 
by a renewal process, but for the sake of our example we shall not do so. 

Let us write a = 1/a, so that a is the mean interarrival time in the 
"thinning Poisson process. " The original process may be viewed as the 
limiting case a = 0+, and as we increase the value of a, the procedure will 
tend to delete clusters selectively, starting with the densest ones first. We 
can therefore expect that the decrease in the tails of the spectral density 
and eventually also for smaller values of w , will reflect (possibly) several 
levels of clustering in the given process. 

For a finite Markov renewal process, it is readily seen that the thinned 
process is itself a Markov renewal process with the same state space. For a 
given positive value of a, the Laplace-Stielt jes transform h(s;a) of the 
transition probability matrix H(x;a) of the thinned Markov renewal process is 
given by 



h(s;a) = [I - h(s + a)]” 1 [h(s) - h (s + a)] 



(18) 



The formula (14) can be numerically implemented for the tranform matrix 
h(s;a) and this for various values of the parameter a. However, the 
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substantial analytic simplifications for the Markovian arrival process make it 
worthwhile to carry out this task for that particular class. This is not only 
because of its prevalence in applications, but also because the physical 
behavior of some MAPs is well understood and they can therefore serve as 
benchmark examples in assessing the utility of the SQSD as a statistical tool. 

For the Markovian arrival process, the matrix h(s;a) is given by 

h(s; a) = a(sl + al- Q) -1 (sI - C) _1 D (19) 

A moment’s reflection shows that h(s;a) is also the transform of the 
transition probability matrix of the MAP with the coefficient matrices and 
of order 2m, given by z: 



Q-al 


al 


, D = 


0 


0 


0 


C 


1 


D 


0 



By elementary calculations, we find that the stationary probability 
vector (a) , © 2 (a)] of the generator Q^ =C^ + D^, is given by 

© (a) = ©D(al - C) \ and © (a) = a©(al - C) 1 (20) 

—l — — 2 — 

The fundamental rate A (a) of the thinned process is therefore given by 

A* (a) = s (a)De = e(I - aC) _1 De (21) 

—2 — 

We note that a graph of the function A (a) is already quite informative, 

particularly if there are intervals over which that function decreases 

* * 

rapidly. The ratio (A ) A (a) is the fraction of transitions of the original 
process that survive in the thinned process. We have found a graph of that 
function useful as an easily computable means of quantifying the "burstiness" 
of an MAP. 

In evaluating the SQSD S(f;a) for increasing values of a, no appreciable 
matrix-analytic simplifications arise from the particular form of the matrices 
C i and D j . Items such as the vector ©(2)D and the matrix (C^-D^) are 

computed once and for all and the further numercial computation consists in 
evaluating the inverses [u 2 I + (C -D J 2 ]” 1 for given values of u> and a. We 

start by computing the SQSD for the given MAP by using formula (17). The 

spectral density should be computed for a sufficiently large value of f, 

so that for larger values of f the function S(f) is negligibly small. To 
determine f is not entirely a hit-or-miss proposition. By using formula 
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(6) to determine the value of f for which the SQSD for a Poisson process of 
rate X decays to a small value, say 10~ 5 , an adequate starting value of f 
is found. In an initial computation, we can divide the interval [0,f ] in, 

max 

say 100 equal parts. Since the spectral density S(f;a) "decays" to zero with 
increasing values of a, that should be adequate to obtain plots for selected 
values of a . These will then suggest other values of a or specific ranges of 
f -values to be explored in further numerical computations. 
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ABSTRACT 

The stress release model is a piecewise linear Markov Process similar to 
the collective risk model, in which the observations consist of the times and 
sizes of the Jumps, which are taken to represent the times and sizes of large 
earthquakes. Earlier studies have shown that the asymptotic behaviour of the 
likelyhood-based tests for this process against a Poisson null hypothesis show 
anomalous behaviour. The present paper develops simulation methods for the 
process and uses them to investigate quantitatively some of the qualitative 
predictions of theoretical studies. In particular it is confirmed that the 
distribution of the likelyhood ratio statistic is non-standard. 
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1. Introduction 

This paper describes two methods for simulating a point process of "stress 
release” or "self-correcting" type, and uses these in an investigation of the distributions of 
parameter estimates and test statistics for this process. The results are of importance in 
fitting the model to data on large-scale historical earthquakes (see, for example, Zheng and 
Vere- Jones (1990)). 

The data to be fitted consists of pairs (t i ,S i ) representing the "occurrence times" 
and "stress releases" from a list of large historical earthquakes. In other words, they can 
be regarded as observations on a marked point process with the times (tj) as points and the 
stress releases (S^ as marks. 

The stress release model is a marked point process with conditional intensity 

function 

(1) 5l(t,S \% t ) = \|/[X(t)] f(S I X(t)) 

where \|f is a monotonic increasing "risk function", X(t) is a piecewise linear Markov 
process used to model the regional stress level, and the probability density f(S lx) describes 
the distribution of the stress release S when the stress level has reached the value x at the 
time when the jump occurs. For the rest of the paper we shall assume that the distributions 
f(S lx) are in fact independent of x, so that the successive values of the stress release form a 
sequence of i.i.d. random variables. 

The process X(t) will be taken to be of the form 

(2) X(t) = X(0) + B[t-CS(t)] 
where B and C are constants and 



1 S { 

i : 0 < t. < t 



S(t) = 
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is the cumulative sum of the stress releases up to time t. To be specific we shall assume 
that the intensity function can be written in the form 

(3) Y(x) = exp{a + bx}. 

Then by combining (1), (2) and (3) we can write the intensity function for the times {tj} in 
the 3-parameter form 

(4) \(i \% t ) = ¥ [X(t)] = exp{ A + B[t - CS(t)] }. 

As for the stress releases, we shall consider three special cases: 

(i) S- constant; 

(ii) S i have a common exponential distribution with density Xe‘** (x ^ 0); 

(iii) S i have a common Pareto distribution with density CqX -8 (x > 1). 

As we shall see, the behaviour of the likelihood ratio statistic is sensitive to the tail 
behaviour of this distribution. 

The case of constant jumps was studied by Isham and Westcott (1979), who 
called the process a "self-correcting” process because of the tendency of the number of 
events N(t) to closely track the target function t/C. Indeed, they established that as t -» ©o, 
VarN(t) remains bounded instead of increasing as 0(t) as in the case of Poisson and 
Renewal processes. This is clear evidence of a high degree of dependence over time, and 
grounds for expecting that the asymptotic inference theory may not apply. This question 
was examined by Ogata and Vere- Jones (1984), who showed that while the asymptotic 
distribution of the Maximum Likelihood (M.L.) parameter estimates was normal, the 
likelihood ratio statistic for testing against the null hypothesis of a simple Poisson process 
was non-standard. Further results are given by Hayashi (1988) and Zheng (1990) 

Inference for the case of variable jumps, while of greater importance in the 
earthquake context, has not been studied theoretically, although Vere-Jones (1988) did 
show that the bounded variance property for N(t) holds only for the case of constant 
jumps. Also, the data sets in practice are usually small (N = 30-100) and the applicability 
of the asymptotic theory is therefore open to doubt 

The present study was undertaken to supplement the theoretical investigations, and 
to give confidence to the application of the model to historical earthquakes. 
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The simulation methods are outlined in Section 2 below; Section 3 treats the 
distribution of the parameter estimates, and Section 4 the distribution of the likelihood ratio 
under the null hypothesis of a simple Poisson process. 



2. Simulation Methods 

We describe two methods, the first a special case of the "inverse transform" 
method, and the second a special case of the "thinning" method, both following the general 
discussion in Ogata (1981). 

The applicability of the inverse transform method to simulating point processes 
rests on the observation that between events the conditional intensity function plays the role 
of the hazard function for the time to the next occurring event (see, for example, Chapter 
13 of Daley and Vere-Jones (1989)). In particular, the time interval £ from an arbitrary 
starting point tg to the next occurring event has a survivor function 

(5) S(x) = P(J; > x I %J = exp[- f X(u 1 tt*(u))du] 

V 

where %*(u) denotes the history of the process up to tg augmented by the assumption that 
no further event occurs between ^ and u, u > t^. 

Substituting the representation (4) for the conditional intensity of the stress release 
model we rind 



(6) J X(u I %*(u))du = J exp{ Ag + Bu}du 

*0 o 

= R 0 (c Bt -1) 

where A 0 = A + B[t 0 - CS(tg)] and Rg = (exp Aq)/B. Note that A 0 and hence R 0 can be 
written down explicitly in terms of the parameters A, B, C and the events (t.,Sj) with 

0 *ti<t 0 . 



The simulation can therefore proceed according to the following algorithm, 
supposing the simulation to terminate after generating a prescribed number of points N or 
after covering a prescribed time interval T. The parameter values A, B, C are supposed 
given. 
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Algorithm 1 

1 . Set t 0 = 0, S(t 0 ) = 0, N(t 0 ) = 0. 

2. Calculate Rg in (6). 

3 . Generate a uniform random variable U, and compute the time £ to the next 
event from 

5 = B' 1 log c [l - Gog e U)/Ro]. 

4. Generate a random variable S from the assumed distribution for the stress 
drops f(x). 

5. Update t 0 , Sftg), N(t 0 ) according to 

t 0 = 

set,,) = s(t 0 ) + s 

N(t 0 ) = N(t 0 ) + 1 

6. Stop if N(tg) > N or t 0 > T. Otherwise record (t 0 ,S) and return to step (2). 

This provides a simple, effective and rapid algorithm for generating successive 
pairs (§.,S.) for interevent times and associated stress drops. 

The thinning approach was first suggested by Lewis and Shedler (1979) in the 
context of non-homogeneous Poisson processes in arbitrary dimensions, where the inverse 
transforms method is inapplicable or involves excessive computations. This is not the case 
here, but nonetheless is may be of interest to examine an alternative method. The 
underlying idea is to first simulate points in a dominating (higher intensity) simple Poisson 
process, and then to successively accept or reject the points in this simulation, with a 
probability of acceptance equal to the ratio of the intensity of the required process to that of 
the dominating process at the point being examined. 

One difficulty with this approach in the present example is that, although the 
conditional intensity (4) is everywhere finite, it has no uniform upper bound. This feature 
complicates the construction of a dominating Poisson process and therefore rather spoils 
the inherent simplicity of the thinning approach. The following version may be suggested. 
It requires an initial choice of the dominating intensity Xq which we discuss below. 
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Algorithm 2 

1. Set t 0 = 0, S(t 0 ) = 0, N(t 0 ) = 0. 

2 . Generate a uniform random variable U and from it an exponential random 
variable by setting 

\ = -(logeUVV 

3. Compute the ratio p = X, (tp + ^V\,, where 

X. (t 0 + £) = exp{A + B[t 0 + £ - CS(t(,)]}. 

(a) If p > 1 evaluate the time for which X.(t Q + ^*) = X Q from 

= {log X 0 -[A + B(t 0 - CS(to))] }/B. 

Set t 0 = t 0 + £*, X^ = 3 X,q. 

Stopift 0 >T. Otherwise return to step 2. 

(b) If p < 1 generate a uniform random variable U. 

(i) If U < p, accept % and generate a random variable S from the 
stress drop distribution. 

If p < x /6 set Xfl = X„ /3, otherwise retain X. Q . 

Update t 0> S(tg), N(t 0 ) according to 

t 0 = l o + ^ 

S(to) = S(to) + S 
N(t 0 ) = N(to) + l 

Stop if N(to) > N or t 0 > T. Otherwise record (t^S) and return 
to step 2. 

(ii) If U > p reject Set t 0 = t 0 + 

Stop if t 0 > T. Otherwise return to step 2. 

A suitable initial value of is about three times the average stress level X, which 
may be found from the following consideration. For large t, assuming the process is in 
equilibrium 
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X(t) = «t-CS(t» = 0(1) 

and so as t — > oo 

S(t)/t C 1 . 

But S(t) » 5 N(t) where "5 is the mean stress drop, and so 
N(t)/t (ScyK 

This means that after the process is started, X(t) will tend to drift upwards (or downwards 
as the case may be) until X(t) is fluctuating around a value in the vicinity of (SC)" 1 . In the 
case of constant jump, only the parameter B determines the local character of the 
realisations. When B is small the risk is close to constant and the process close to Poisson; 
when B is large, the intervals are close to regular, the risk starting at a low value just after 
an event and reaching larger and larger values as the interval increases. 

In the case of variable jump size an important role is played also by the tails of the 
jump distribution. If the distribution is highly peaked around a mean value, the behaviour 
will be similar to the constant jump case. If occasional very large jumps can occur, the 
process will exhibit a quiescent period after the large jump, followed by more frequent 
smaller events until the next large jump occurs. Some examples of trajectories for different 
combinations of parameter values and distributions are shown in Figure l(a)-(d). 



3. Distribution of Parameter Estimates 

In this section the simulations are used to examine the properties of the parameter 
estimates obtained by standard likelihood maximization procedures. The three main 
purposes are: 

(i) to check that the algorithms being used are working correcdy; 

(ii) to confirm the theoretical results, in particular consistency and the asymptotic 
forms of variances, obtained for the maximum likelihood estimates by Ogata and 
Vere- Jones (1984); 

(iii) to check the extent of small sample deviations from the asymptotic forms. 

The log likelihood here has the characteristic point process form 

log L = £ log X(t.) - f X.(u)du 

0<t<T 1 o 



( 7 ) 
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where X is given by (4) and contains the three unknown parameters A, B, C. The integral 
can be evaluated explicidy as a series of exponentials over the intervals between events. It 
can therefore be used diiecdy as the target function in a standard optimization procedure. 
In our case the NAG routine (E04JAF) was used, which requires no derivatives to be 
given. Quadruple precision was used in all evaluations of the target function, to avoid 
rounding errors in these calculations causing instability in the iterations and hence 
producing error messages. Note that with this particular form of parameterization, a 
unique solution is guaranteed by convexity (cf. Proposition 2 of Ogata and Vere-Jones 
(1984)). 



Since this paper treats only the constant jump case, this was chosen for the 
simulation, with unit jump size and parameter values A = 3, B = 2 and C = 1. The 
parameterization in (4) above differs slighdy from that used in the earlier paper, viz 

(8) X(x) = exp {a + |3[t - p 0 N(t)] + ytyT} 



where p 0 = C is the true value. The asymptotic relationships between the standard errors 
of the two sets of parameter estimates are given by 



(9) 



o, = 

A a 



a* = a A 

B P 

a A = r 1 . 

C 1 



the standard errors on the right being given by 



Var[T 1/2 (6t, 



J 



-1 



where 



J = 



Aj V 
V w |v 

, JLtj lv JLtt 
\^2 U 2 V 3 u y 



and 
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i T 

U = lim ^ f A.(u)du 

t-^-t 

, T 

V = lim f X(u)A.(u)du 
T-> ” T o 
T 

W = lim i r X 2 (u)X(u)du 
T-*~ T o 

with ^(u) given by (8) (or equivalently (4)) and X(t) = t - p 0 N(t). 

To use these results it is first necessary to carry out a preliminary simulation to 
evaluate U, V, W. This was done in two ways, firstly by directly tracking the trajectory of 
X(t) over a long period (T = 10,000) and substituting into the above limits, and secondly 
by using the expressions given in Vere- Jones and Ogata (1984) and Ogata and Vere- Jones 
(1984) for U, V and W as expectations with respect to a discrete skeleton of the process 
X(t). This required finding the stationary distribution of the discrete skeleton and 
evaluating the expectations of the necessary functionals. The two approaches gave values 
which agreed to within 1 % which was considered adequate. Averaging the two sets of 
values gives 



ful 

V 




f 1.003^ 
-1.267 


; J = 


f 1.00 -1.27 .50 ^ 
-1.27 1.81 -.635 


; J'* = 


f 11.5 6.0 -5.8^ 
6.0 4.8 0.1 ! 


IwJ 




1.809 J 




^ .50 -.635 .334 J 




■ 

U\ 

oo 

O 

VO 



As a further check, we know that N(t)/t — » y 1 = 1 from the bounded variance 
property, while E(U) = E[X(t)] = E[N(t)/t] in the stationary regime. Thus in theory U = 1. 
Also, from the form of J, the (2,3) position of J' 1 must be 0. 

Substituting finally from (9) we obtain for the limit standard errors and correl- 
ations, the predictions 



T 1/2 fr A -4 3.4 -4 0.8 

T 1/2 ft B -> 2.2 fr AC -4 0.5 

T 3/2 fr c -4 1.73 fr BC -4 0 

The process was then simulated 100 times for T = 32, 64,1 2 8,256, 512, para- 
meter values being estimated for each run. The means, standard errors and correlations are 
set out in the table below. 
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TABLE 1 

Simulation results and theoretical limits for parameter estimates 





32 


Length of Simulation 
64 128 


256 


512 


Limit 


Means 


A 


3.30 


3.19 


3.11 


3.02 


3.04 


3.0 


B 


2.21 


2.14 


2.08 


2.02 


2.02 


2.0 


C 


1.001 


1.0002 


1.0003 


0.99997 


1.00003 


1.0 


Standard Errors 


T 1/2 S a 


2.60 


3.57 


3.41 


3.78 


3.89 


3.4 


T 1/2 Sd 


1.78 


2.23 


1.90 


2.47 


2.44 


2.2 


t 3/2 s® 


1.59 


1.84 


1.98 


1.91 


1.62 


1.73 


Correlations 


r AB 


.67 


.76 


.68 


.82 


.88 


0.8 


r AC 


.49 


.53 


.67 


.51 


.58 


0.5 


r BC 


-.18 


-.07 


.01 


-.00 


.20 


0.0 



The results overall are reassuring. The procedures are producing sensible results, 
broadly consistent with the theory. There is evidence of positive bias in the estimates of A 
and B for small T, but it is at worst of the order of 10% and decreases with T. The results 
for the asymptotic variances are borne out well by the simulations, even for quite modest 
sample sizes. 



4. Distribution of the Likelihood Ratio Statistic 

The discussion in Section 3 of Ogata and Vere-Jones (1984) shows that in the case 
of constant jumps, the likelihood ratio statistic 

(10) A = 21og(t 1 /t 0 ) 

where Lj is the likelihood (7) for the stress release model using the parameters estimated 
by maximum likelihood, and L 0 is the likelihood for the simple Poisson model with 
estimated rate, has a non-standard asymptotic distribution. Instead of following the % 2 
distribution on 2-degrees of freedom, it has the distribution of a certain quadratic functional 
of standard Brownian motion on [0,1]. This phenomenon arises from the long-term 
memory property of the stress release model with constant jumps, which makes possible 
the "bounded variance" property of the trajectories N(t), and causes the Hessian matrix 
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(2nd derivative of the log likelihood) to converge to a random not a constant limit. The 
case of variable jump size has not been investigated theoretically, but it seems likely that 
similar phenomena will appear also. 

Since the form of the distribution is crucial to the tests and selection procedures 
(AIC) used with the historical earthquake data, it is important to gain an idea of the extent 
of the departures from the standard % 2 form. Since, moreover, the true distribution is 
either unknown (variable jump case) or hard to evaluate (constant jump case), and is in any 
case unknown if the sample sizes are small, simulation offers the only practicable route 
towards ascertaining the size of these departures. In this section we outline the results of a 
preliminary investigation of this kind, for various combinations of sample size and jump 
distribution. 

On general grounds, we know that the distribution of the statistic (10) should be 
independent of the scales of measurement used for the time and stress dimensions. It 
must, therefore, be independent of the size of the jump in the case of constant jumps, and 
of the parameter X in the case the jumps follow exponential distribution 1 - e’** . The case 
of most practical interest, however, is that of a Pareto distribution with power law tails, 
say 



( 11 ) 



l-F(x) 



C 0 x e 

.0 



x£x 0 
x< x Q 



since a Pareto distribution for stresses corresponds to the exponential distribution 
("Gutenberg-Richter relation") generally used for earthquake magnitudes. Here the ratio 
may be influenced by the shape parameter 0, but not by the threshold value x Q . 



In all a total of 8 simulation runs were performed for constant jumps, exponential 
jumps, Pareto jumps with 0 = 2 (corresponding to b = 1.5), 0 = 1 (corresponding to b = 
0.75 approximately in the Gutenberg Richter relation), each run for T = 32 and T = 64 time 
units. Note that in the second case the sufficient conditions for ergodicity in Zheng (1990) 
are not satisfied since the jump distribution has infinite mean. 



Approximate values of some summary statistics for the distributions of A are 
shown in Table 2, and some sample density functions shown in Figure 2. They are based 
on 500 simulations each for T = 64, T = 32. Note that in the estimation part of the 
program, the parameters B and C of the stress release model were constrained to be 
non-negative. 
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TABLE 2 

Simulation results for the distribution of the likelihood ratio statistic A 





Constant Jump 
T = 32 T = 64 


Exponential Jump 
T = 32 T = 64 


Pareto Jump (1) 
T = 32 T = 64 


Pareto Jump (2) 
T = 32 T = 64 




Mean 


3.7210 


3.7862 


2.5139 


2.4854 


2.1669 


2.2539 


1.1385 


1.4079 


2.000 


Median 


3.5633 


3.4666 


2.3644 


2.0579 


■Rrmn 


1.8566 


0.7686 


0.9583 


1.386 


LQ 


2.3446 


2.1500 


1.2299 


0.9783 






0.2315 


0.3443 


0.575 


UQ 


4.9671 


5.0424 


3.4754 


3.3832 




3.1736 


1.7059 


2.0545 


2.770 


90%ile 


6.2100 


6.7771 


4.5638 


5.1747 


4.7719 


4.8911 


2.7216 


3.4282 


4.605 



From these results a number of surprising if tentative conclusions can be drawn. 

(a) In all cases the use of the % 2 results would be misleading. The exceedance 
probabilities (p-values) of a given observed value are higher than those for the % 2 
form for constant and exponential jumps and for the first Pareto form, but lower 
than the % 2 value for the second Pareto form 

(b) The discrepancy between the real and % 2 form is at its most acute for the case of 
constant jumps. 

(c) The discrepancies relate both to the shape and the scale of the distribution, and are 
suggestive of an "equivalent degrees of freedom" which decreases as the tail of the 
jump distribution becomes thicker. 

It would clearly be highly desirable to have some theoretical underpinning of these 

results, and an attempt to provide such undeipinning is the next stage in our project. 



Acknowledgements 

The authors are grateful for the opportunity to contribute to this celebration 
volume, and in so doing to recognize the outstanding contributions Professor S K 
Srinivasan has made to stochastic modelling in many fields, above all in point processes 
and their applications. 

The work was commenced while the first author was a Visiting Fellow at Victoria 
University on sabbatical leave from the University of Malaya. The support of both 
institutions is gratefully acknowledged. 






23 



List of Figures 

Figure 1 Simulated trajectories of stress release models. 

a. Constant Jumps 

b. Exponential Jumps 

Figure 2 Histograms for the distribution of the likelihood ratio, fitting 
a stress release model to Poisson data. 

a. Constant Jumps 

b. Pareto Jumps 0 =1.5 

All graphs show values of A along the x-axis and frequency up the y-axis. 
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POSITIVE DEFINITE FUNCTIONS IN QUANTUM MECHANICS AND IN TURBULENCE 



J. BASS 

Universite Pierre et Marie Curie, Paris 
ABSTRACT 

Positive-definite functions have various representations which can be 
interpreted in probability theory, in quantam mechanics, and in turbulence 
theory. The equivalence of these representations is studied. The nature of 
the characteristic functions in quantum mechanics is discussed. In 
turbulence, time -correlations are positive-definite functions, but 
space -correlations have an ambiguous structure. In the annex, the properties 
of pseudo-random functions, which are convenient representations of 
turbulence, sire summarized. 

INTRODUCTION 

In probability theory, a probability space is given. Probability measures 
over this probability space are chosen and measurable functions (random 
variables) are defined. 

In the applications to physics, it happens that "probability laws" are 
introduced, by means of their Fourier transforms (characteristic functions), 
which are continuous complex- valued functions of one or several real 
variables. The characteristic functions may have various forms, each of them 
being associated to a particular sort of "statistical mechanics". But the 
characteristic functions which are employed in a given statistical mechanics 
may not be compatible with an unique underlying probability space. Quantum 
mechanics offers an example of this situation. 

The aim of this paper is to give a review of the main representations of 
characteristic functions, to examine their connections and to discuss the 
question of compatiability. 

1 - Positive-Definite and Characteristic Functions 

We consider complex- valued functions f(A) which are continuous and 
positive definite. For ary choice of integeral n, complex C »C 2 » ...,C n and 
real X ,A , . . . , A they are defined by the condition : 

12 n 
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Z C C f (A - A ) > 0 (1) 

k 1 k 1 - 

These functions have several representations and interpretations as 
characteristic functions and as correlation functions. 

I - According to the fundamental theorem of Bochner f is the Fourier transform 
of a positive bounded measure m : 



00 

f(A) = S exp(iAs) dm(s) (2) 

“00 

f(A)/f(0) is the characteristic function of some probability law. 

II - Let H be a Hilbert space and 0 an element of H such that II 011 = 1. Let A 
be a hermitian operator over H, bounded or sometimes not bounded. If < > 
denotes the scalar product in H, the expression 

f(A) = < exp(iAA)0 ,0 > (3) 

is a positive-definite function. If it is continuous, it is the 

characteristic function of a probability law. 

In quantum mechanics, this probability law is interpreted as the law of 
the quantity associated with the operator A , in the state 0, by a "rule of 
correspondance " . 

III - Let a be a real -valued function defined over IR, or over IR + . The mean 

values 

T T 

f (A) = lim 2 ^ f ex P (iAa(t)) dt, or lim i J exp (iAa(t)) dt (4) 

T— >oo ^ -T T — >00 o 

are positive-definite functions, if they exist. 

They are used in the deterministic theory of turbulence and can be 

interpreted as characteristic functions of probability laws. 

IV - Let Q be a probability space and s an element of Q. We chose a 

probability measure m(s). Let p(s,t) be a stationary random function of 

second order. Then the mathematical expectation 

r 00 - 

E p(s,t) p(s,t+x) = p(s,t) p(s,t+x) dm(s) (5) 

-00 

is a positive-definite function of x (independent of t). It is interpreted as 
the covariance of the random function p. This function is the Fourier 

transform of a "spectral measure". 
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V - Let g be a complex valued function such that, for every x , 

T 

f(x) = lim F g(t) g(t+x) dt (6) 

T— >oo -T 

exists. Then f(x) is a positive definite function. It is interpreted as the 
correlation function of the function g. 

The proofs of all these results are almost obvious. We see that there 
exist two main interpretations of a positive definite function : 

i) characteristic function of a probability law ; 

ii) covariance or correlation of a function, random or not random. 

2. TWO-DIMENSIONAL CHARACTERISTIC FUNCTIONS 

The definitions in section 1 are for the case of one-dimensional 
characteristic functions. They are easily extended to several dimensions; 
but, in the cases II and III (quantum mechanics and turbulence), this extension 
needs precautions. 

a) Let A,B be two hermitian operators and consider 

f(A,p) = < exp(i(AA+pB) 0, 0 > (7) 

If A and B commute, it is easy to verify that f(A,p) is positive-definite. It 

is the Fourier transform of a probability measure m(s,s').If A and B do not 

commute, this result need not be true; f(A,p) is not generally a 

characteristic function. Only some particular points 0 in the Hilbert space 
have this property [11, 12]. 

2 

A classical example can be given when H is a L -space, and 
A = Multiplication by x (operator of position) 

B = y ”j x (operator of momentum) 

It is found that 



f ( A, #i) = J exp(iAs) ^ [ s ~^] dS 



(8) 



f is the Fourier transform of 



kj exp(_i,ls) *[ s + *§] * ( s - ^ 



(9) 
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This is Wigner’s pseudo-density function (called Wigner-Bass pseudo-density by 
L. de Broglie in [8]). 

It is easy to find functions such that (9) is not everywhere positive. 
For instance, this situation is verified when 

0(s) = -J— for | s | <1 

= 0 for |s| >1 

This result shows that the probabilistic structure of quantum mechanics 
is not compatible with a unique underlying probability space. This remark 
appears more clearly in terms of operators. The set of all hermit i an 
operators is a vector space. Any family of commutative operators belongs to a 
subspace. But two non-commutative operators belong to independent subspaces 
and the intersection of these subspaces is itself a smaller subspace of the 
whole space. 

2 

Let us take the example of the hamiltonian operator H in a L -space. It 
does not commute with x or with p = j- The family of unitary operators 

exp(iHt) (an abelian group) is used for generating the solution of Schrodinger 
equation. We have 

0(x,t) = exp (iHt) 0(x,O) (10) 

<exp(iHt)0,0 > can be identified to be the characteristic function of the 
hamiltonian operator, the time t playing the role of the variable X . 

In the space of operators, the trajectory of exp(iHt) belongs to a 
subspace e of commutative operators. But this subspace is different from the 
spaces c 2 » of x or of p. The situation is roughly represented by fig. 1 
e i ,C 2 ,€ 3 have onl y the common point 1 (operator identity). 




Cl 

^2 



Figure 



x 
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b) The characteristic functions of pairs of functions give rise to very 
similar structure. The characteristic function of the pair of functions 
a(t),0(t) is the mean value 



f(A,p) = lim f exp i[Aa(t) + p0(t)]dt 

T-*» -T 



( 11 ) 



if this limit exists. For instance, the characteristic function of 
(cos t, sin t) is equal to 



lim 

T-*» 



J exp i[Acos(t) + 0sin(t)]dt = Ja 2 + p‘ 
-t ° 



( 12 ) 



Of course, the associated probability is concentrated over the unit circle. 
But it happens that the mean value (11) does not exist, even when a(t) and 
0(t) possesses good characteristic functions. For instance, the formula 

i fexpUA log t) dt = (13) 

O 

shows that log t has no characteristic function. But it is easy to see that 
t+logt has a characteristic function. That is due to the fact that t/logt is 
an increasing and invertible function, not essentially different from t. As t 
possesses a (discontinuous) characteristic function, At+p(t+log t) has 
characteristic function if A+p * 0, but not if p = - A. 



We say that two function a, 0 (real or complex) are comparable if the 
mean value of a0 exists. Comparability of functions is analogous to 
commutativity of operators. A family of comparable functions constitutes a 
vector space, a sub-set of the set E of all functions possessing a quadratic 
mean value. The set E is not a vector space. But it is itself a subset of a 
vector space, namely the Marcinkiewiz space , which is the space of complex 
valued functions a such that 



1 r 

lim sup == |a(t ) | dt < co 
T— >oo -T 



(14) 



E can be dissociated into vector sub-spaces of comparable functions. As an 
example, the functions exp(it) and exp[i(t+log t)] have quadratic mean values, 
but are not comparable. 



3. CONSTRUCTION OF FUNCTIONS WITH A GIVEN CHARACTERISTIC FUNCTION 

We restrict ourselves to the case of a characteristic function which is 
the Fourier transform of an absolutely continuous measure. We are given a 
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characteristic function f . According to Bochner’s theorem, we have 

r 00 

f(A) = exptiAs <p(s)]ds , <p £ 0 (15) 

-00 

We can write 

f ( A ) = <exp i\sJ<p(s) , 4 tp(s)> (15) 

where ^ € L 2 . Therefore there is certainly a representation of f under the 
form (II) of section 1. 

Is it possible to represent f as 



1 f T iAa(t) 
lim = e dt 

T-*» o 



(we use 




rather than 




for practical reasons). 



(17) 



For this purpose, we shall introduce step functions, defined by uniformly 
distributed sequences (see annex). 



Let z be a uniformly distributed sequence over ]0, 1[. If A is a 

n 

function defined over ]0,1[ and Riemann-integrable, we know that (Weyl’s 
theorem) . 



Therefore 



1 N r 1 

lim ± £ A(z n ) = J A(z) dz 
n-x» n=o n o 



lim 

N oo 



1 

N 



N 

E 



n=o 



exp(iAA(z n )) = J exp(iAA(z)) dz 

o 



(18) 



(19) 



Let a(t) be a function equal to A(z ) for n<t<n+l). Then left-side member of 

n 

(19) is exactly 



lim 
T oo 




iAa(t) .. 
dt 



The question reduces to finding A such that 

f (A) = J l e UA(z) dz (20) 

O 

It is well known that it is possible. For instance we choose for A a function 
increasing from -oo to +oo when 0<z<l, and introduce the inverse function 



such that 
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A(z) = u, z = A j (u) 

If A is differentiable, 

1 

00 

f < A) = | e 1Au A'(u) du (21) 

where A' is the Fourier transform of f. 

This is nothing but a constructive interpretation of Bochner’s theorem. 

Let us now study the case of a two-dimensional characteristic function 
c*(A, p). 

First Method: We introduce a two dimensional sequence (x ,y ) uniformly 

n n 

distributed in the square (0,1) x (0,1). We suppose that the two sequences 
(x ) and (y ) are independent (see annex). 

n n 

In order to solve the equation 

i r T 

f ( A, p) = lim ^ J ex P i(Aa(t)+fi0(t))dt (22) 

T co o 

we represent a and |3 by step functions, namely 

a(t) = A(x ,y ), 0(t) = B(x ,y ), for n<t<n+l (23) 

n n n n 

We have 

1 N 

f (A, p) = lim rr V exp i[AA(x ,y ) + jliB(x ,y )] 

N L n n n n- 

T oo n = o 



=JJ e l(AA(x,y) + MB(x,y)) dxdy (24) 

c 

where C is the square (0,l)x(0,l). (24) has the ordinary form of a 

characteristic function. If X and Y are two independent random variables, 
uniformly distributed in the square C , f is the characteristic function of 
the pair of random variables A(X,Y), B(X,Y) . 

If the mapping 

5 = A(x,y) 7) = B(x,y) 

is invertible, and if J is the jacobian, we have 

m.n) = [ f J_ d?d7l (25) 

J J IJI 

and -1— is the probability density corresponding to f. 

PI 
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We have only to choose A and B such that 



D( A, B) 
D(x, y) ' 



( 26 ) 



Q A 

If A is a function of x only and not of y , then J = ~ — 5 — and 

n n n OX dy 

exp [iAA(x )] is the characterist ic function of the function a. This procedure 
n 

gives an inductive method for constructing an n-dimensional characteristic 



function. 



Second method - We shall show that the introduction of two sequences x ,y is 

1 1 ■ n n 

not necessary. Let z be a unique uniformly distributed sequence, and 

n 

a(t) = A(z ), 0(t) = B(z ) for n < t <n+l (27) 

n n 

Can we find A,B such that 

f(X, „) = 1 im I + * lB(x » )) (28) 

N oo n = 0 



- J 1 e 1(XA(x ’ y) + ** B(x,y)) dxdy ? (29) 

O 

Here A(X) and B(X) are two random variables, functions of X, uniformly 
distributed over (0,1). They are not independent. 

Let X be a random variable, uniformly distributed over ( 0 , 1 ). Let us 
assume that it is possible to find two functions U(X), V(X) such that the two 
random variables U(X), V(X) are independent and uniformly distributed over 
(0,1). We define A(x) and B(x) by 



Then 



A(x) = F i [U(x),V(x)], B(x) = F 2 [U(x),V(x)] 
f(A,H) = J exp i[XF i (U(x), V(x)] + pFJUfx) , V(x)] dx 

O 



(30) 



(31) 



f(A,p) is the stochastic mean value (expectation) of the random variable 
AF i (U, V) + #iF 2 (U, V), defined over the probability space constituted by the 
interval (0,1). Therefore 



f(A.„) = JJ e 1 (^ 1 (u.v)+HF a (u.v)) du dy (32) 

c 

where C is the square (0,1) x (0,1) . 

As in the first method, we see that it is possible to choose F and F g in 
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such a way that f is a given function. 

We must now prove that is possible to find two functions U,V 
over (0,1), such that the random variables U(x), V(x) are independent, 
let us recall the definition of Rademacher functions. We represent 
number x, 0 ^ x ^ 1 by its infinite binary expansion : 



00 




a 

k 



a = 0 or 1 

k 



Instead of (33), we prefer writing 



oo r 

1 - 2x = V — — , with r = 1 - 2a , r = 1 or -1 
L 0 k k k k 

k= 1 2 

r is a function of x, called Rademacher function, 
k 

It is easy to verify that r k (x) is a step function, such that 

r (x) = (-I)”, -Ell- < x < , p =1,2 2k 

k 2 k 2 k 

One has the formula [9]: 

I exp i £ [A^r^tx)] dx = ^ J exp [iA^Cx)] <*x 
o ^ k =i k k J o 

We give the proof in the case where n=2, which corresponds 

two-dimensional characteristic functions. We have to compare 

f = J exp [i(Ar j (x) + pr 2 (x))] dx 

o 

f^= J exp [iAr i (x)]dx ; f g = J exp [ifir 2 (x)]dx. 

o o 

It is easy to verify that 

f = j |exp[i (A+jii)] + exp[i(A-/n)] + exp[-i(A-p)] + exp[-i (A+p)]j 
f i = 2 t ex P^^^ + ©xp(-iA)] = cos A 



f = - [exp(ip) + exp(-ijit) + exp(ip) +exp(-ip)] = cos p. 

2 4 



defined 
First 
a real 

(33) 

(34) 

(35) 

(36) 
to the 
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Therefore 







1 

1 

1 

1 

1 


r i 


jl/2 

1 

1 

! 


jl X 

i 

i 

i 


i r i 

i i : r * 

1 1 


%\ ’«>! %! 

1 ‘ 

i ! L 


r x 

i 



exp 



[i(Ar i (x) 



|xr 2 (x))] 



dx 



-r 



exp 



[iAr^x)] 



dx 



exp 



[ i M^r ^ ( x ) ] dx 



(37) 



This formula signifies that the random variables r^Cx), r 2 (x) are independent. 

The functions r j , r g have the values -1, 1 . But it is clear that the property 

(37) is still valid for the functions a , a , which take the values 0, 1 . 

1 2 



The functions a^(X), a (X), resulting from the binary expansion of X, 
constitute two independant random variables. They can be used as functions 
U(X), V(X) . 



4 - SOME QUESTIONS CONCERNING CORRELATION FUNCTIONS. 

We have seen that positive-definite functions can be interpreted not only 
as characteristic functions of a probability law but also as correlation 
functions. 

We shall show that, in quantum mechanics as well as in turbulence, this 
notion of correlation leads to some contradictions. 

a) In quantum mechanics , a hermit ian operator A, associated with a state <p t 
generates a probability law and plays the role of a random variable. As $ is 
a function of time, this random variable seems to be a random function of 
time. At two times t, t' , we have the pairs (A, 0(t)) and (A,^(t')), 
equivalent to a pair of random variables X(t), X(t'). Is it possible to 
attribute a significance to the mean value of the product X(t).X(t')?. 

As the mean values are defined by a unique 0, we must reduce (A,0(t')) to 
the form (B, tfr(t)), in such a way that (A, \jj( t')) and (B, ^r(t)) correspond to 
the same statistics ; they must have the same probability law, i.e. the same 
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characteristic function : 

<exp( iAA)0(t' ) , 0(t')> = <exp( iAB)0(t ) , 0(t)> (38) 

But 0 satisfies Schrodinger equation. If H is the hamiltonian operator, the 
evolution of 0 is given by 

0(t') = exp[iH(t'-t)0(t)] (39) 

Therefore 

<exp( iAA)0(t' ) , 0(t' ) )> = <exp(iAA) exp[iH(t'-t)0(t)], exp[iH(t'-t )0(t)]> (40) 

As exp(iAH) is a unitary operator, (40) is verified if 

exp[-iH(t'-t)] exp(iAA) exp[iH(t' -t )] = exp(iAB) (41) 



This equality for all values of A is equivalent to 

B = exp[-iH(t'-t)]. A. exp[iH(t'-t )] (42) 

B is well defined. But in general this operator B does not commute with A. 
Thus even if we interpret (A, 0(t)) as a random function, we cannot attribute 
any significance to the covariance of this function. 

This remark makes it clear that quantum mechanics is different from 
classical statistical mechanics. The probability laws in quantum mechanics 
cannot be related to a unique probability space. 

b) Let us now consider the problem of turbulence . The experiments in fully 
developed turbulence strongly suggest that turbulence must be represented by 
deterministic functions, namely pseudo-random functions of time (J.BASS 
[3,4]). The idea of a deterministic representation of turbulence is also 
contained in the new theories of dynamic systems, which concern the birth 
and evolution of turbulence, rather than the study of a stationary state of 
flow [7] . 

The notion of time correlation of turbulence suffers no special 
difficulty. Let us introduce a symbol for the representation of time mean 
values. We represent by Mf the mean value 



Mf = lim 

T-*co 




(43) 



If u(x,t) is the velocity of the fluid, as a function of time at a given point 
x, the time correlation is given by 
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F(x,x) = M u(x, t )u(x, t+x) dt (44) 

But the "turbulencists" introduced, around the years 1930-40, the notion of 
space correlation, between u(x,t) and u(x+h,t) , measured at two 
different points at the same time. It is still a time mean value, namely 



R(x,h,t) = M u(x,t) u(x+h,t) dt (45) 

In general, it is a function of x and h . It has no property of homogeneity, 
i.e. stationarity in space. As the fluid is confined in a bounded box, 
clearly its properties at the center of the box are different from those near 
the walls. 



Consequently, there is no reason for R(x,h, t) to be the Fourier transform 
of a spectral measure. In particular R is not an even function of h. 

Of course, we improve the structure by transforming the correlation 
function (bilinear covariance) into a correlation coefficient : 



r(x, h, t ) 



R(x, h, t ) 



J Mu 2 (x,t) J Mu 2 (x+h,t) 



(46) 



The expansion of r in the vicinity of h = 0 is 



r(x, h, t ) 




Mu 2 Mu ' 2 - ( Muu ' ) 2 1 



(Mu 2 ) 2 



(47) 



By Schwarz’s inequality, we verify that the coefficient of h 2 is positive. 
The experimental evidence shows that: 



1. The time corrleation is rapidly decreasing (with small oscillations) and 
tends to zero when x -» ». In other words, the velocity is a pseudo-random 
function (terminology introduced by J.BASS ([2, 3, 4, 5 ]). 

2. The space correlation has, when h increases, a very similar shape. It 
approaches zero rapidly and this fact is not related to the presence of the 
bounder ies. 



A function u(x,t) can then represent the velocity of a turbulent fluid if 
it satisfies the following requirements : 

i) At every point it is a pseudo-random function of time ; 

ii) Its space correlation - a time mean value - is rapidly decreasing and 
approaches zero when h increases. 

We say that such a function is a turbulent function . 
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It is easy to give examples of turbulent functions. The most elementary one 
is a function defined by 



u(x,t) = f(x+at) (48) 

where f is a pseudo-random function and a a constant. If y is the correlation 
function of f, the time correlation of u is y(ax). The space correlation is 
y(h). This function is pseudo-random in spacee and in time . 

More generally, let f be a sequence of independent pseudo-random 

n 

functions : 



M f (t) f (t+x) = 0 for n * p , 

n p 




= y (x) for n = p 


(49) 



(correlation function of f ). 

n 

The existencee of independent functions is explained in [5] . 

If c n (x) is a sequence of functions of space, the function 

u(x,t) = Z c (x) f (x+at) (50) 

n n 

is a turbulent function (for an infinite series, trivial conditions of 
convergence are necessary). In fact, the mean value of u(x, t )u(x, t+x) and 
u(x, t)u(x+h, t) are 

M u(x, t )u(x, t+x) = Z c (x)c (x)M f (x+at)f (x+at+ax) 

n p n p 

=Z|c(x)| 2 y(ax) (51) 

1 n 1 n 

M u(x, t)u(x+h, t ) = Z c (x)c (x+h)Mf (x+at)f (x+h+at) 

n p n p 

= Z c (x)c (x+h) y (h) (52) 

n n n 

The decay of the space correlation when h increases is related to the 
structure of f , and does not depend on the coefficients c (x). But in this 

n n 

case there is no homogeneity in space. 

As a conclusion any theoretical approach of turbulence must take into 
account the fact that, for a stabilized turbulence, the functions which 
represent the velocity of the flow must be not only pseudo-random in time, but 
"turbulent J ‘ in space and time. 
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Annex - Uniformly distributed sequences. 

A sequence of real numbers is uniformly distributed over (0,1) if the 
following properties are satisfied : 

i) 0 < z < 1 . 
n 

Let N' be the number of points z ,z , ....z such that a<z <b, where (a,b) is 

1 2 n n 

an arbitrary sub-inteerval of (0,1). Then 

. . N' 

lim = b-a . 

v N 
N*oo 

An analogous definition is available for a sequence of points belonging to the 
cube (0, l) p of IR P . (a,b) is replaced by an arbitrary paral lelopiped. 

A sequence of numbers is p-uniformly distributed if the sequence 

{z z > is uniformly distributed in the cube (0,l) p . A theorem of 
n+p-l , n+p 

H. weyl (1916-see [5]) gives a condition for a sequence to be uniformly 

distributed :a sequence z is uniformly distributed if, for any Riemann- 

n 

integrable function <f > , one has 



lim 

N-*» 




(54) 



An equivalent condition is that, 

N 



lim - 

N 

N-)oo 



z 

1=1 



for any integer l different from zero, 
exp(2imz ) = 0 

n 



(55) 



Let P(z) = a z p + a z p_1 +...+ a be a polynomial of degree p £ l with real 

0 1 p 

coefficients. If a is an irrational number, the sequence {p('n), modulo 1}, 
o 

is p-uniformly distributed . 

The sequence {© n , modulo 1} is uniformly distributed for almost every © . 
But we don't know explicitly the values of © having this property. Only 

counter examples are known: for © = — the sequence © n is not uniformly 

distributed modulo 1. 

If the sequence z is 2-uniformly distributed, the function equal to 

n 

0 (z ) for n<t<n+l, and 0 for t<0, is pseudo-random. Its correlation function 
n 

is equal to A(l-|x|) U(x) where 

U(x) = 1 if | x | ^ 1 , and 0 if |x| £ 1 
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A = J 1 0 ( z ) | 2 dz. 

o 

Two uniformly distributed sequences z , s are independent if the sequence 

n n 

lz^ + l's^ is uniformly distributed, whatever be the integers i,i' not 

simultaneously equal to zero. It is equivalent to saying that 

1 N 

lim rr V exp [i ( iz +I's )] = 0. 

„ N 4-1 n n 

N->oo n= 1 

If z^ and s^ are independent, 

1 N 

lim - Y 0 (z ,s ) = || 0 (z,s)dz ds 

^ n = 1 n n J 

where C is the square ( 0 , 1 ) 2 . 

Examples of independent sequences : two polynomials P,Q of different degrees. 

Two polynomials p = a Q z p +...,Q = b Q z p + ...of the same degree, when the 

irrational numbers a , b are independent (i.e. Aa +ub =y, A,u,r integers, 

o o o o 

implies A = p = r = 0; V2 and i/3 are independent but V2 and l+\/2 are not). 
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ABSTRACT 

The quantum Input-Output formalism of Gardiner and Collett is used to derive the rate equation describing the 
photoelectron counting statistics of an electromagnetic field output from a high-Q optical cavity. The evolution of the 
cavity field is assumed to be governed by an arbitrary Hamiltonian for a single mode of the field. The resulting 
equation is shown to conform with one previously derived on intuitive grounds, using purely population statistical 
arguments. As an explicit example, the solution for a freely-evolving cavity field is computed. 



Glossary of Principal Terms Used in Text 



a(tUZ in (tla 0 Jt)) 

d, D 

FKE 

f (*; t) 

h(s'; t) 

H 

h 

10 

MGF 

m 

n 

m 

N 

P(n; t) 

P(n, N; t ) 

Qte t) 

QJfi 

Q (*,*'; t) 
y 

K /*, v 

v 

v 

P 

0) 



cavity, {input, output} field operator. 

detection operators. 

forward Kolmogorov equation. 

(l-jy*>. 

(1— IfS 

Hamiltonian operator. 

Planck’s constant. 

Input-Output. 

moment generating function. 

number of photons (individuals) in cavity population at time /. 
cavity photon number operator. 

number of photoelectron counts registered by detector in interval [0, t). 

output number operator. 

probability distribution for n. 

joint probability distribution for n and N. 

MGF of p(n ; t), with expansion parameter s. 
initial cavity MGF. 

MGF of P(n, N; t), with expansion parameters s and s'. 
cavity loss rate. 

population birth, death, and immigration rates, respectively. 

detector quantum efficiency. 

vacuum field for detector. 

density operator. 

cavity mode angular frequency. 
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1. INTRODUCTION 

Many quantum optical systems can be described using classical population statistics: indeed, Srinivasan has 
extensively employed the formalism of classical point processes to analyse a range of systems representing cavity 
radiation models 1 . When an electromagnetic field is describable in terms of photon number, the evolution of the 
photon statistics can often conform to those of a population model, and visualisation of the behaviour of the system at 
the microscopic level is greatly simplified. 

For such systems it is comparatively straightforward to include a mechanism for direct detection of the field 
within the population model 2 ; one such process has been dubbed population monitoring 3 . To date, however, the 
justification for this detection formulation has involved purely intuitive arguments, although for some particular 
examples of field statistics the population monitoring approach has yielded results equivalent to those obtained by 
purely quantum field theoretic methods 2 . It would be preferable to demonstrate the equivalence ab initio. 

In the following we shall use the quantum Input-Output (10) methods formulated by Gardiner and Collett 4, 5 to 
derive an equation for the photoelectron statistics for a detected field from a high-Q cavity, and show that, in the 
situations where the cavity field alone can be described as a simple population process, the photoelectron statistics are 
equivalent to those obtained by population monitoring. 



2. POPULATION MONITORING 

Consider a population containing n(t) individuals at time t, where n(t) is a Markov process. Let p{n\ t) represent 
the probability of there being n individuals in the population at time t. Now define the moment generating function 
(MGF) for this distribution as Q(s\ t), where 

00 

t) = £ Pto 0 (i-*)" • (2- 1 ) 

n=0 

We shall assume that the evolution of p(n; t) is governed by a forward Kolmogorov equation (FKE) which can be 
expressed in terms of the MGF as follows, 

gf ( *° = <P(S. d/ds)Q(s-.t). (2 2) 

Here, ip(s, d/ds) is a differential operator appropriate to a given system; for example, if the population evolution is 
governed by homogeneous birth, death, and immigration processes, with rates A, (jl> and v, respectively, the FKE is 

p(n ; t) = A(n-l) p(n- 1; t) + //(«+!) P(n+U t) + vp(n-\\ t) - (v + (A +p)n) p{n\ t) (2.3) 

and from (2.1) and (2.2) the corresponding function tp is 

ip(s, d/ds) = — vs — ps + ?is(l—s) ^ . (2.4) 

For a single mode of a quantised cavity electric field, n{t) represents the corresponding photon number. 
Detection of the field leaking out of the cavity may be incorporated into this picture by extending p{n\ t) to the joint 
distribution P{n y N\ t), defined as the joint probability that there are n photons in the cavity at time t, and that N(t) 
photoelectrons have been counted at the detector in the semi-open interval [0, f), ( i.e. excluding any count at the point 
t). The manner in which N(t) is determined from the process for n(t) is defined by population monitoring 3 . In this, the 
photons are assumed to leak from the cavity in accordance with a death process of rate y, say, and are registered at 
the detector with probability (quantum efficiency) rj. The associated extended FKE for P(n,N ; t) contains terms, 
associated with the cavity processes, similar to those in the equation for p(n; t), but now merely appended by the 
parameter N y as these processes involve no change in the photoelectron number. The detection processes above 
introduce the additional contributions to the FKE for P(n, N\ t ): 6 

,< rt ’ N '>*) = Vy (»+l) P(*+l, JV-1; t) + (n+1) P(n+ 1, N; t) - yn P{n y N\ t ). (2.5) 

aaaitlonai 

Upon introduction of the MGF Q for the joint distribution P(n, N\ t), 
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QHs.s'-.t) = £ £ P(n.N\ t) (1-s)" (1-s'f , (2.6) 

n=0 N=0 

the general FKE for P(n, N\ t) is transformed to 

|p = tp(s, a/a s) q - y s %* + ws'lp • (2.7) 

Setting s'= 0 in Q (s t s'\ t) recovers the MGF Q(s\ t) for the cavity population alone. Comparison of eqs. (2.2) and (2.7) 
clearly reveals that the effect of population monitoring on the undetected cavity population is merely to introduce an 
additional death process of rate y. The term ijys'dQ/ds in (2.7) is responsible for the specific detection characteristics 2 . 

In practice, it has often been found that if (2.2) has an exact analytic solution, then so also does (2.7). For a 
Markov process, the solution to (2.7) provides all the information needed to determine any multilinear moment of the 
counting distribution. The population monitoring procedure has also been applied by Srinivasan to several 
non-Markovian population systems 1 . 



3. THE IO FORMALISM 

Gardiner and Collett have developed a formalism describing a single mode of a high-Q cavity quantum electric 

4 5 

field coupled to the continuum of field modes outside the cavity ’ . In the following we consider only a one-sided 
cavity, in which one mirror is transmitting and the other totally reflecting. For an isolated system the field mode is 
described by the Heisenberg picture annihilation and creation operators a(t) and afy), which evolve according to the 
Heisenberg equation 

kt) = - i/h [ 2(f), H 1 . (3.1) 

where H is the given system Hamiltonian, and h is Planck’s constant divided by 2 n. When coupled to the outside 
world, however, the mode evolves according to the modified equation 

3(0 = - i/h l a(t), H ] — y/2 2(0 + . (3.2) 

Here a ln (t) is the Fourier transform of an input external field mode operator, and y is the cavity damping constant, 
proportional to the transmission coefficient of the input-output mirror. The operators a in (t), possess the 

commutator 

|2 to (0, 2 JV)| = a(r-t') . (3.3) 

In addition, a further external operator represents the output from the cavity. This operator has a commutator 

similar to that for ajj), 

l 2 «W. 2<J(01 = Mf-f) , (3.4) 

and is related to a in (t) and d(t) through the boundary condition 

2(0 = y ' 1/2 ( a ln (i) + Z 0 Jt ) ) . (3.5) 

Causality is assured by demanding that the following commutators hold: 

[ F,(2t(0, 2(0), F/ajO, 2^(0 ] = 0 , ?>t, (3.6) 

[ F 3 (2t(0, 2(0), f/2 0u ,(r'), 2 ol >') ] = 0 . t>? , (3.7) 

(where F v F v F y and jF 4 are given functions), implying that an input field at a particular time can have no effect 
upon the cavity field at an earlier time, and that a cavity field can have no effect upon an output field at an earlier 
time. Equations (3.2), (3.5), (3.6) and (3.7) ensure that unitarity is maintained in this scheme, i.e. 

[2(0,2^01 = 1 v t. (3.8) 

More general system operators c(() (e.g. functions of 3(0, 2^(0) evolve according to the modified Heisenberg 
equation 
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c(0 = — «/* [ c«, 5 1 - { |c(t), 2 t (i)| (y/2 2(r) - y 1/2 2,„(0 ) 

- (y/2 Sty) - 7 1/2 Sj(0 ) [2(0 . 2(01 [ , (3.9) 

and the operators ( y/2 a{t) — y 1/2 a in (t)) and (y/2 a(t) — y l/2 o out (t)) can be shown to commute with all system 
operators at time t. 5 



4. PHOTODETECTION IN THE IO FORMALISM 

In order to relate the IO formalism sketched above to the photodetection problem described earlier, we shall 
assume that the cavity in question is one-sided and that no light is input to the cavity. Thus, the input state is assumed 
to consist solely of the vacuum. The internal cavity field evolves in accordance with eq. (3.2), and the output field a out 
interacts with the detector, which registers the photoelectron counts. A general detector quantum efficiency rj is 
included by noting that the relevant operator giving the statistics of the counting process is d(t), given by the unitary 
transformation 

2(0 = TJ 1 ' 2 a o Jt) + (1-rjf 2 v(t) , (4.1) 

where v(r) is a vacuum field, contributing nothing to normally-ordered moments of d{t) and d^(t), and commuting with 
all other operators 7 . Moments of the photocounting distribution for the interval [0, t) are given by moments of the 
detection operator D(t), where 

m = JV 5V') 5(,') , ( 4.2) 

so that, for example, the mean-square of the number of photoelectrons N(t) counted in the interval [0, t) is 

m,)f) = Tr[p(6(0) 2 ]. (4.3) 

where Tr denotes the quantum mechanical trace over the subsequent operators, and p is the relevant density operator. 

We need to construct from these operators an MGF that generates such moments. The obvious choice 
corresponding to eq. (2.6) is 

o («,*';<) = ((l-sy^a-s'y* 0 ), (4.4) 

where 

2(0 = ’ aV) 2(0 (4.5) 

is the photon number operator for the cavity field. 

In principle, Q (s, s'\ t) in (4.4) generates the joint distribution P(n, N\ t). This quantity is not uniquely defined 
quantum mechanically, however, until it can be demonstrated that the operators n(t) and D(t) commute. Now, since 
v(/') commutes with all other operators, and D(t) consists of operators a ouX {t') and v(f') with t' < t, it follows from 
the causality relations in (3.7) that D(t) commutes with all system (cavity) operators at time t. Hence 

[n(0,£(*)l = 0, (4.6) 

and the MGF in (4.4) is uniquely defined. Note that (4.6) implies that the number of photons in the cavity at time t 
and the number of photoelectrons counted up to time t are simultaneously measurable. 



5. EFFECT OF QUANTUM EFFICIENCY 

In order to derive an evolution equation for Q(s t s'\ t) we must eliminate the operator 3(t) in favour of the 
operator a ouJ (t). To do this, we shall employ the normal-order theorem for continuous-mode operators 8 . For arbitrary 
function q(t), this reads 

exp ( j d /' q(t ') 3V) 2(0 ) = : exp ( j dt' ( exp q{t') - 1 ) d\t’) 3(*')) : , 



(5.1) 




